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1  Introduction 


The  development  of  future  aerospace  vehicles  requires  the  optimization  and  coupling  of 
diverse  technologies  in  order  to  operate  under  stringent  conditions  near  the  performance 
envelope.  The  design  of  such  vehicles  is  critically  dependent  on  simulation  tools  which 
integrate  several  fundamental  disciplines  and  which  necessarily  exploit  the  power  of 
advanced  parallel  computing  platforms.  The  motivation  to  develop  high-fidelity  algorithms 
derives  from  considerations  in  various  areas  of  current  interest.  For  example,  in  fluid 
dynamics,  the  effort  to  minimize  the  impact  of  aerodynamically  generated  sound  on 
communities  and  structures  has  necessitated  large-scale  time-domain  aeroacoustics 
computations.1  Another  field  where  numerically  intensive  wave-propagation  phenomena 
dominate  is  in  electromagnetic  wave  scattering,  which  is  crucial  in  determination  of  radar 
signature  and  in  biomedical  modeling.  Similarly,  in  the  field  of  turbulence  simulation, 
methods  to  overcome  the  traditional  limitations  of  Reynolds-averaged  approaches  have 
focused  on  direct-  and  large-eddy  simulation  techniques  which  seek  to  resolve  the 
multispectral,  highly  nonlinear  stochastic  phenomena  from  computationally  challenging 
first  principles.  Endeavors  which  combine  fluid  dynamics  with  other  disciplines  also  yield  a 
large  and  typically  stiff  equation  set  whose  numerical  solution  mandates  the  development 
and  utilization  of  high-resolution  techniques.  One  example  is  the  complex 
flow-acoustic-structural  interaction  which  gives  rise  to  practically  significant  phenomena 
including  panel  flutter,  transition  delay,  drag  and  noise  reduction,  turbulent  flow  sensing 
and  modification,  and  control  of  unsteady  separation  by  means  of  deforming  flexible 
surfaces.  Another  area  of  great  current  focus,  magnetogasdynamic  flow  control,  derives 
from  the  integration  of  fluid  dynamics  with  electromagnetics  and  models  of  thermochemical 
nonequilibrium.  The  potential  benefits  include  drastic  reductions  of  thermomechanical 
fluid  load  and  control  of  transition,  turbulence  and  separation,  without  the  deployment  of 
moving  surfaces,  which  is  a  very  attractive  option  at  hypersonic  speeds. 

The  equations  governing  the  behavior  of  the  above  phenomena  of  interest  exhibit  a  wide 
range  of  mathematical  properties  and  are  described  in  Section  2.1.  The  high-resolution 
method  required  to  resolve  the  stiff  spatiotemporal  scales  must  also,  therefore,  exhibit 
highly  flexible  traits  necessary  for  uniform  applicability,  with  little  regard  to  the  structure 
of  the  equation  set.  Additional  difficulties  arise  from  the  fact  that  modern  computational 
methods  require  focus  on  practical  problems  with  their  attendant  geometrical  complexity. 
Configurations  of  interest  are  typically  three-dimensional  and  numerical  simulation  is 
greatly  facilitated  by  curvilinear  or  body-fitted  coordinate  systems.  Since  the  phenomena 
are  typically  unsteady,  and  boundaries  may  be  flexible,  the  ability  to  treat  deforming 
meshes  is  also  essential. 

A  survey  of  the  literature  indicates  a  clear  requirement  for  new  versatile  high-resolution 
methodologies  to  address  the  sophisticated  design  space  described  above.  In  this  report,  a 
powerful  combination  of  very  high-order  finite-difference  schemes  is  developed,  which  fulfils 
this  need.  The  method  is  based  on  compact  (Pade-type)  formulas2,3  of  fourth-  and 
sixth-order  accuracy,  which  are  employed  to  form  the  various  derivatives  (Section  2.2.1). 
Although  these  schemes  have  a  long  history,  their  use  has  been  restricted  to  geometrically 
simple  situations  in  which  uniform  meshes  suffice  and  boundaries  may  be  treated 
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simplistically.  Application  to  practical  multidisciplinary  problems  has  languished  due  to 
their  extreme  susceptibility  to  numerical  instabilities  which  are  initiated  at  mesh 
nonuniformities  and  at  truncated  boundaries  and  which  subsequently  grow  unbounded.  To 
overcome  these  difficulties,  the  present  effort  supplements  the  difference  schemes  with 
filters.3,4  This  component  is  also  Pade-type  as  described  in  Section  2.2.2  and  is  utilized  in  a 
post-processing  framework  to  remove  high-frequency  oscillations  which  are  not  supportable 
on  the  meshes  employed. 

Several  adequate  time-marching  methods  are  available  for  use  with  spatially  high-order 
accurate  methods.  As  described  in  Section  2.2.3,  here  we  adopt  the  classical  fourth-order 
four-stage  Runge-Kutta  method  for  primary  use  in  wave  propagation  dominated  situations. 
On  fine  and  stretched  meshes  of  the  type  employed  in  wall-bounded  flows,  this  explicit 
scheme  is  unacceptably  restricted  in  allowable  time-step-size.  In  such  cases,  an  implicit,  up 
to  third-order  accurate  approximately  factored  scheme  is  employed  within  a  sub-iteration 
framework. 

The  difficulties  associated  with  the  use  of  higher  order  schemes  on  general  curvilinear 
meshes  are  not  restricted  to  the  appearance  of  numerical  instabilities.  A  potentially 
debilitating  source  of  error  is  associated  with  the  failure  to  preserve  freestream  when 
metrics  are  evaluated  in  the  standard  manner.  In  Section  2.2.4,  elegant  procedures  are 
developed  to  eliminate  this  source  of  inaccuracy  without  significant  computational  penalty. 

The  rapid  gain  in  popularity  of  relatively  cheap  parallel  computing  platforms  effectively 
demands  that  new  computational  techniques  be  capable  of  exploiting  distributed 
processing  strategies.  Pade-type  methods  are  not  trivially  adaptable  to  such  approaches 
due  to  their  implicit  nature.  In  Section  2.2.5,  a  domain-decomposition  technique  is  derived 
for  arbitrary  meshes.  The  key  element  is  a  multipoint  overlap  procedure  to  treat  the 
artifical  boundaries  created  at  domain  interfaces  in  general  curvilinear  grids. 

To  validate  the  method  and  to  highlight  its  properties,  Section  2.3  discusses  several 
applications  in  some  of  the  areas  of  interest  noted  above.  For  the  purposes  of  brevity,  these 
represent  only  a  small  subset  of  the  problems  examined  but  nonetheless  emphasize  the 
versatility  of  the  approach. 

In  addition  to  broad  spatio-temporal  scales,  high-fidelity  methods  are  also  required  for 
multifluid  problems  associated  with  thermochemical  nonequilibrium.  For  this  class  of 
problems,  an  effort  is  described  to  understand  the  interaction  between  various  molecular 
kinetic  modes  such  as  vibration- vibration  (V-V)and  vibration-translation  (V-T),  which 
have  a  profound  impact  on  the  performance  of  high-speed  devices.  Vibrational  freezing  in 
the  nozzle  of  a  high-speed  propulsion  engine  represents  wasted  energy  that  is  not  recovered 
as  thrust  work.  A  prerequisite  to  minimizing  this  loss  is  to  develop  an  understanding  of  the 
dominant  thermochemical  processes.  In  Section  3,  advances  in  Master  and  extended 
Navier-Stokes  simulation  techniques,  defining  the  physics  of  the  vibrational  manifold  are 
documented. 
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2  High-Order  Methodology  for  Fluid  Dynamics  and 
Electromagnetics 

2.1  Governing  Equations 

In  order  to  highlight  the  multidisciplinary  nature  of  the  compact-difference-based 
technique,  consider  a  generic  governing  equation  system  written  in  curvilinear  coordinates 
(f  =  V,  z,  t),  rj  =  7}(x,  y,  z,t)  C  =  CO,  y,z,t),r  =  t )  in  nondimensional, 
strong-conservation  form: 


d_  /Q\  dF_  dG  dH  _S 
dr  V  J )  dp  d(  J 


(1) 


Here  Q  denotes  the  solution  vector,  J  —  d  (£,  p,  C,  r)  /d  (x,  y,  z,  t )  is  the  transformation 
Jacobian  and  S  is  a  source  vector.  The  form  of  the  flux  vectors  F,  G  and  H  depends  on 
the  nature  of  the  problem  being  solved  and  may,  in  general  be  written  in  terms  of  the 
corresponding  Cartesian  flux  vectors,  F,  G  and  H .  For  example, 


F  =  -j{ixF  +  iyG  +  izH). 


(2) 


Various  physical  models  may  be  obtained  by  suitably  choosing  the  vectors  in  Equation  1. 
For  brevity,  in  the  following  we  focus  primarily  on  the  vectors  Q ,  F  and  F  with  the 
understanding  that  other  vectors  may  be  written  in  an  analogous  fashion  or  obtained  from 
the  citations  provided.  The  Euler  equations  are  obtained  (see  e.g.,  Ref.  5  for  details)  when 
Q  is  the  five-component  vector,  {p,  pu,  pv,  pw,  pE},  S  =  0  and  F  —  Fj  where 


Fj 


pu 

pu 2  +  p 
puv 
puw 

(. pE  +p)u 


(3) 


For  subsequent  reference,  F 
by 


Fj  for  the  Euler  equations  in  curvilinear  coordinates  is  given 


Fj 


pU  a 

pull  +  CxP 

PVU  +  iyP 

pwu + izP  ^ 
(pE  +  p)U  -  ip 


(4) 


In  these  equations,  p  is  the  density,  u,v,  w  are  Cartesian  components  of  the  velocity,  p  is 
the  static  pressure,  E  is  the  total  energy,  p/ (7  —  1  )p  +  ( u 2  +  v 2  +  w2)/ 2,  U  is  the 
contravariant  velocity  vector,  U  —  i  +  £xu  +  £yv  +  £,zw,  and  i  =  =  J^1^,  ■  ■  ■ 

denote  the  metrics  of  the  coordinate  transformation. 
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The  Navier-Stokes  equations  may  be  obtained  by  adding  components  associated  with  the 
viscous  fluxes,  i.e.,  F  =  Fj  +  Fy  where 


Fv  =  { 


-rT 


Re  TxV 

J-r 
Re  lrr 


'  xz 


(' UTXX  +  VTxy  +  WTXZ)  +  7— 


OTe  i UT- 


~8T 


(7— l)PrM2Re  dx 


(5) 


Here,  Txx,rxy, . . .  denote  components  of  the  stress  tensor  (see  Ref.  5  for  detailed  formulas), 
Re  is  the  Reynolds  number,  M  and  Pr  are  the  Mach  and  Prandtl  numbers,  T  is  the 
temperature  and  7  is  the  ratio  of  specific  heats. 

To  express  Maxwell’s  equations  in  the  template  of  Equation  1,  Q  may  be  set  to, 

Q  {Exi  Fy:  EZ1  Bx:  By ,  Bz }  and 


F 


< 


0 


0 


Ez_ 


>  • 


l^m 


(6) 


where  /im  and  e  are  the  permeability  and  permittivity  of  free  space  and  E—{Ex,Ey,Ez} 
and  B={BX,  By,  Bz }  are  the  electric  field  and  magnetic  induction  vectors  respectively. 


The  governing  equations  of  magnetogasdynamics  are  obtained  by  combining  Maxwell’s  and 
Navier-Stokes  equations,  together  with  Ohm’s  law,  under  the  assumption  that  relativistic 
effects  may  be  ignored  and  that  the  displacement  current  is  small.6  In  this  case,  Q  in 
Equation  1  is  set  to  the  eight-component  vector  { p ,  pu,  pv,  pw,  pE ,  Bx ,  By,  Bz}  where  the 
total  energy  E  now  includes  components  associated  with  the  magnetic  field, 

E  —  p/ ((7  —  l)p)  +  ( u 2  +  v2  +  w2)/ 2  +  RbB  ■  B / (2 p,mp).  Correspondingly, 


F 


pu 

pu2  +  P-  Rb^Bx 
puv  -  Rb^Bx 
puw  —  Rbj^Bx 
(pE  +  P)u-  R„  (U  ■  £)  B, 
0 

uBy  —  vBx 
uBz  —  wBx 


(7) 
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where  P  =  p  +  is  the  sum  of  the  static  ( p )  and  magnetic  pressures,  and 

z/irn.  v  ' 


Fv  = 


1  1  (d_By_  _  d_  BaA 

Rea  a  \dxfim  dy  fim  ) 

i  i  _  jlmA 

Rea  a  \dx  fim  dz  fim  ) 


(8) 


Here  a  is  the  electrical  conductivity,  Rh  =  - jyp1 -  is  the  magnetic  pressure  number 

PrefUiefy,mref  K 

and  Ra  =  LUreffimrefo-ref  is  the  magnetic  Reynolds  number. 

It  is  important  to  recognize  that  in  writing  the  governing  equations  in  strong  conservation 
form,  Equation  1,  the  following  metric  identities  have  been  invoked: 


h  =  {ix)t  +  (Vx)v  +  (Cz)c  =  0-  (9) 

h  =  {iy)e,  +  ( Vy )rl  +  (4)c  =  0-  (10) 

h  =  (4)?  +  (Vz)v  +  (4)c  =  0-  (11) 

I4  =  (1/ J)r  +  (4)?  +  +  (Ct)c  =  O'  (12) 


The  importance  of  satisfying  these  identities  in  the  present  higher  order  finite-difference 
numerical  scheme  is  detailed  in  Section  2.2.4. 


2.2  Numerical  Technique 

As  noted  earlier,  a  finite-difference  approach  is  adopted  to  solve  the  discrete  equivalent  of 
the  governing  equations.  Thus,  the  values  of  the  solution  vector  are  localized  in  a  pointwise 
sense  at  each  node  of  the  mesh.  This  considerably  simplifies  the  procedure  of  formal 
extension  to  higher  order  accuracy  and  contrasts  with  the  finite-volume  approach  where 
cell  averaged  quantities  are  typically  employed. 


2.2.1  Difference  Scheme 


The  approach  is  illustrated  by  considering  a  generic  scalar  pointwise  discrete  quantity  0. 

At  different  stages  of  the  solution  procedure,  this  could  be  a  metric,  a  primitive  variable  or 
a  component  of  a  flux  vector.  In  all  cases,  the  required  spatial  derivative  (f)'  is  obtained  in 
any  transformed  coordinate  direction  by  solving  the  tridiagonal  system: 


,  it  ,  r  ,/  Ai+ 2  _  4>i-2  4>i+ 1  _  4>i- 1 

r  (pi-1  +  </>*  +  T0i+1  =  b - - - +  a - - - . 


(13) 
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where  T,  a  and  b  determine  the  spatial  properties  of  the  algorithm.  Variations  of  this 
formula  have  been  described  in  several  references  (see  for  example  Ref.  3).  It  encompasses 
a  family  of  schemes  ranging  in  accuracy  from  the  standard  three-point,  second  order 
explicit  method  (E2),  with  T  =  0,  a  —  1,  b  —  0  to  the  compact  five  point,  sixth  order 
algorithm  (C6)  with  T  =  -I  a  —  b  —  5.  Another  scheme  of  interest  is  the  three-point, 
fourth  order  method  C4  with  T  =  a  =  ^,  b  =  0. 

The  scheme  of  Equation  13  cannot  be  applied  at  a  few  points  near  each  boundary  where 
the  stencil  protrudes  the  domain.  Here,  special  one  sided  boundary  schemes,  such  as  those 
described  in  Refs.  7  and  8  are  employed.  A  compendium  is  provided  in  Ref.  9. 

To  compute  the  residual,  the  derivatives  of  the  inviscid  fluxes  are  obtained  by  first  forming 
the  fluxes  at  the  nodes  and  subsequently  differentiating  each  component  with  the  above 
formulas.  Both  viscous  and  nonideal  magnetogasdynamic  terms,  if  present,  are  obtained  by 
first  computing  derivatives  of  the  primitive  variables.  The  components  of  the  viscous  flux 
are  then  constructed  at  each  node  and  differentiated  by  a  second  application  of  the  same 
scheme.  Although  this  approach  is  not  as  accurate  as  that  in  which  a  Pade-type  scheme  is 
employed  directly  for  the  second-derivative,  it  is  significantly  cheaper  to  implement  in 
curvilinear  coordinates.  A  previous  effort10  demonstrates  that  successive  differentiation 
yields  a  stable  method  in  conjunction  with  the  added  filtering  component  which  is 
described  next. 

2.2.2  Filtering  Scheme 

Compact-difference  schemes  are  susceptible  to  numerical  instability  due  to  their  centered, 
nondiffusive  nature.  Most  problems  of  practical  interest  exhibit  nonlinear  behavior,  require 
domain  truncation  and  nonperiodic  boundary  conditions.  Additionally,  the  regions  of 
interest  are  usually  discretized  with  nonuniform  meshes.  Each  of  these  complications  can 
introduce  spurious  oscillations  which  are  not  damped  by  centered  schemes.  Thus,  in 
practical  situations,  an  additional  component  is  required  in  the  method  to  enforce 
numerical  stability.  The  approach  taken  here  is  to  eliminate  spurious  frequencies  through  a 
low-pass  filtering  procedure. 

Denoting  a  typical  component  of  the  solution  vector  by  (f),  filtered  values  0  satisfy, 

1  +  +  Olf(j)i+ 1  =  E^=0—  {(j)i+n  +  4>i-n)  ■  (14) 

This  template  has  been  obtained  from  Refs.  3  and4.  Equation  14  has  N  +  2  unknowns, 
af,ao,ai,...a,N,  which  can  be  derived  from  Fourier  and  Taylor-series  analyses.  The 
constraints  employed  include  the  consideration  that  the  highest  (or  odd-even)  mode  be 
completely  annihilated  while  maintaining  higher  order  accuracy.  The  coefficients 
a0,  ai, . . .  ajv  are  expressed  in  terms  of  atf  which  is  considered  to  be  a  free  parameter  but 
must  lie  in  the  interval  —0.5  <  af  <  0.5.  Within  this  range,  af  controls  the  strength  of  the 
filter  and  as  it  is  reduced,  a  wider  band  of  frequencies  at  the  higher  end  of  the  spectrum  is 
preferentially  damped.  Computations  on  a  range  of  2-D  and  3-D  problems  suggest  that  on 
meshes  of  reasonable  quality,  a  value  0.3  <  ay  <  0.5  is  appropriate.  Only  in  cases  where 
the  mesh  is  of  extremely  poor  quality,  if  it  contains  metric  discontinuities  for  example,  a 
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lower  value  of  aj  ~  0.1  may  be  required.  The  impact  of  filtering  on  the  fully  discretized 
1-D  advection  equation  with  periodic  end  conditions  has  been  examined  in  Ref.  11. 

The  relatively  large  stencil  of  very  high  order  filters  requires  special  formulations  at  several 
points  near  the  boundaries.  For  instance,  the  tenth  order  interior  filter  requires  an  eleven 
point  stencil  and  thus  cannot  be  applied  at  the  near-boundary  points  1, . . . ,  5  and 
correspondingly  at  IL  —  4, ... ,  IL  where  it  protrudes  the  boundary.  The  dependent 
variables  at  points  1  and  IL  are  specified  explicitly  through  the  boundary  conditions  and 
are  not  filtered.  At  the  remaining  near-boundary  points,  two  approaches  have  been 
developed.  In  the  first,  the  order  of  accuracy  is  reduced  on  approaching  the  boundary  to 
the  extent  necessary  for  a  lower  order  centered  formula  to  be  applicable.  Since  meshes  are 
usually  finer  near  the  boundary,  no  significant  accuracy  degradation  is  observed 
particularly  if  a  corresponding  adjustment  (or  optimization)  of  the  value  of  a?f  is  applied.10 
In  the  second  method,  introduced  in  Ref.  12  higher  order  one  sided  formulas  have  been 
developed.  Although  the  spectral  characteristics  of  these  approaches  are  nonideal,  several 
applications  have  demonstrated  the  elegance  of  these  techniques.13 

After  each  time-step,  the  filter  is  applied  sequentially  to  the  updated  solution  vector,  Q ,  in 
each  of  the  three  spatial  directions.  Many  other  strategies  may  be  suggested  but  have  not 
yet  been  explored  in  the  present  effort. 


2.2.3  Time-integration  Schemes 


The  spatial  differencing  scheme  has  been  coupled  with  both  explicit  and  implicit 
time-integration  schemes.  The  classical  fourth  order  four-stage  Runge-Kutta  method 
(RK4),  implemented  in  low-storage  form,  is  utilized  primarily  for  wave  propagation 
problems.  Details  may  be  found  in  Ref.  14  and  are  not  repeated  here. 

Explicit  methods  are  not  suitable  however  when  very  fine  and  highly-stretched  meshes  are 
employed  as  in  wall  bounded  flows.  To  address  this  deficiency,  a  robust  and  highly  accurate 
implicit  procedure  has  been  implemented  through  the  formula: 


J~1P+1 

J-lP+1 

J~1P+1 


+  ^ArSf  ($y 

+  <m^2)  (W) 

+  <M  r8f  (19 


JP+1X 

JP+1X 

A  Q 


=  —d>1  At 


J 


-1P+1  (l-j»QP-(l+2</>)Q"+0Q" 


+Qp[1/J)tp  +  5^Fp) 
Sv  (< Gp )  +  (HP 


At 


+ 


-l 


(15) 


where  OF jdQ  etc.  are  flux  Jacobians,  5  represents  the  spatial  difference  operator  and 
A Q  =  Qp+1  —  Qp.  The  method  combines  the  approximate  factorization  procedure  of 
Ref.  15  with  the  diagonalized  simplification  of  Ref.  16  and,  to  reduce  errors  associated  with 
these  simplifications,  a  sub-iteration  strategy  is  employed.  Thus,  for  the  first  subiteration, 
p  =  1,  Qp  =  Qn  and  as  p  — >  oo,  Qp  — >  Qn+1.  Typically,  three  subiterations  are  applied  per 
time  step.  Note  that  while  the  derivatives  of  the  flux  Jacobians  have  been  obtained  to 
second-order  accuracy  (denoted  with  the  superscript  (2)),  those  on  the  right  hand  side,  i.e., 
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in  the  residual,  are  evaluated  with  the  compact- difference  based  higher  order  method.  In 
the  case  of  the  magnetogasdynamic  equations,  a  loosely  coupled  approach  which  separates 
the  flow  and  magnetic  held  vectors  has  been  described  in  Ref.  17.  Nonlinear  artificial 
dissipation  terms,18  not  explicitly  shown  in  Equation  15,  are  appended  to  the  implicit 
operator  to  enhance  stability.  Although  higher  order  accuracy  in  time  is  achievable  with 
relative  ease  at  the  expense  of  increased  storage  requirements,  a  range  of  numerical 
experiments  suggests  that  second-order  accuracy  in  time  is  adequate  for  the  problems 
considered.19 

2.2.4  Spatial  and  Temporal  Metric  Evaluation 

As  noted  earlier,  several  metric  identities  (Equations  9  through  12)  have  been  invoked  to 
obtain  the  strong  conservation  form  of  the  governing  equations,  Equation  1.  In  a 
finite-difference  formulation,  spatial  and  temporal  metrics  must  be  evaluated  carefully  in  a 
manner  designed  to  satisfy  these  identities  numerically  so  that  freestream  preservation 
errors  are  controlled  to  within  the  machine  accuracy  range. 

For  the  standard  second-order  explicit  method  (E2),  the  straightforward  averaging 
procedure  described  in  Ref.  20  is  adequate.  However,  this  technique  is  not  appropriate  for 
higher  order  methods  -  it  has  been  shown  in  Ref.  21  that  freestream  preservation  errors 
accumulate  to  the  point  where  the  advantages  of  the  new  methodology  are  lost.  A  similar 
loss  of  accuracy  occurs  when  analytic  metrics  are  employed  (if  available). 

To  address  this  important  issue,  we  utilize  the  formulation  developed  in  the  context  of 
lower  order  methods22  in  which  the  metric  relation,  for  example, 

L  =  yvZ(  ~  Vc zv  (16) 

is  evaluated  by  considering  the  equivalent  “conservative”  form: 

( yriZ)(  (y(Z)r]  (17) 

As  demonstrated  in  Ref.21,  this  straightforward  modification  assures  that  the  identities  I\ 
through  J3  (Equation  9  through  Equation  11)  are  satisfied. 

To  satisfy  the  identity  J4,  Equation  12,  also  known  as  the  Geometric  Conservation  Law,  to 
numerical  accuracy,  the  remarkably  simple  approach19  utilizes  1 4  to  establish  the  value  for 
(1  /J)ti  he., 

(1/ J)T  =  —[(&)£  +  (fjt)v  +  (Ct)d-  (18) 

where 

[•£t(i-L:)  T  yriHy)  T  ^r(Cz)] 

fit  =  -[xT(fjX)  +  VriVy)  +  Zr(fjz)\  (19) 

C t  [^(Cr)  T  y-riCy )  +  ^r(Cz)]- 

This  strategy  compensates  for  errors  in  evaluating  time  metrics  and  grid  speeds  which  can 
then  be  obtained  either  by  analytic  means,  if  known,  or  high  order  discrete 
approximations. 19 


Swirl  velocity 


(a)  Vorticity  contours  at  three  instants  with  C4-  (b)  Swirl  velocity  at  T  =  12  along  Y  =  0  with 

i^lO0'4  scheme  several  schemes 

Figure  1:  Single-Domain  Calculation  on  30  x  60  Mesh 

2.2.5  Interface  Treatment  in  Multidomain  Calculations 

As  noted  earlier,  the  major  issue  in  multidomain  computations  is  the  treatment  of 
interfaces  between  domains.  In  the  present  work,  communication  between  adjacent  meshes 
is  conducted  through  finite-size  overlaps.  To  explore  the  issues  involved,  consider  the 
unsteady  inviscid  flow  due  to  a  convecting  vortex  in  an  otherwise  uniform  flow  at  a 
freestream  Mach  number  M ^  =  0.1.  The  initial  condition  is  determined  from  the 
theoretical  solution  given  in  Ref.  23.  The  nondimensional  vortex  strength  parameter, 
C/(UooR)  is  set  to  0.02  where  R  is  the  nominal  vortex  radius.  The  domain  considered 
extends  —6<X  =  x/R<  18,  —6<Y  =  y/R<6.  The  baseline  mesh  consists  of  60  x  30 
uniformly  spaced  points,  with  a  spacing  AX  =  AY'  =  0.41,  which  is  deliberately  coarse 
enough  to  distinguish  between  the  accuracy  of  the  various  schemes.  For  this  problem,  the 
RK4  scheme  is  employed  exclusively,  with  At  =  0.005  which  is  sufficiently  small  to 
guarantee  time-step-size  independent  results  on  each  of  the  meshes  considered. 

As  a  reference,  the  C4  compact  scheme  coupled  with  an  interior  filter  of  tenth  order 
accuracy,  RIO0  4  is  chosen  for  analysis.  Vorticity  contours  obtained  in  a  single-domain 
calculation  with  periodic  end  conditions  are  shown  in  Figure  1(a)  at  three  instants, 

T  =  0,  6  and  12.  At  time  t  =  0,  the  vortex  center  is  located  at  X  =  0  and,  with  increasing 
time,  it  convects  to  the  right  at  a  nondimensional  velocity  of  unity.  The  superiority  of  C4 
over  the  explicit  E4  scheme  is  clearly  demonstrated  in  Figure  1(b)  which  shows  the  swirl 
velocity  along  Y  =  0  at  T  =  12.  Results  with  two  lower  order  schemes  are  also  plotted. 

The  centered  E2  scheme  shows  primarily  dispersive  error,  whereas  the  popular  third-order 
upwind-biased  Roe  scheme24,25  is  highly  diffusive. 
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Mesh  1  5— p,t  py^rl^p  Mesh  2  , 


^  i  i  i  i  i  i 

i  i  i  i  i  i 


IL— 5  IL-4  IL-3  IL-2  IL-1  IL  \ 

•  Meshl  O  Mesh  2  Arrows  show 
w  .  information  transfer 

Figure  2:  Schematic  of  Mesh  Overlap  with  Five  Points 


2.2.6  Mesh  Overlap  Configuration 

In  the  first  set  of  multidomain  calculations,  a  vertical  interface  is  introduced  into  the  above 
mesh  as  shown  schematically  in  Figure  2.  The  computation  is  started  with  the  vortex 
located  at  the  center  of  Mesh  1  (solid  lines).  With  increasing  time,  the  vortex  convects  into 
Mesh  2  (dashed  lines)  after  passing  through  the  interface.  The  meshes  communicate 
through  an  overlap  region.  Although  the  overlap  points  are  collocated,  they  have  been 
shown  slightly  staggered  in  Figure  2  for  clarity. 

At  every  time-step,  the  solution  is  advanced  independently  in  each  domain  with  individual 
interior  and  boundary  formulas  in  the  same  manner  as  in  single-domain  computations. 

Data  is  exchanged  between  adjacent  domains  at  the  end  of  each  sub-iteration  of  the 
implicit  scheme  (or  each  stage  of  RK4),  as  well  as  after  each  application  of  the  filter.  To 
illustrate  data  transfer  between  the  domains,  details  of  a  five-point  vertical  overlap  is  also 
shown  in  Figure  2.  Each  vertical  line  is  denoted  by  its  i-index.  The  values  at  points  one 
and  two  of  Mesh  2  are  set  to  be  identically  equal  to  the  corresponding  updated  values  at 
points  IL  —  4  and  IL  —  3  of  Mesh  1.  Similarly,  reciprocal  information  is  transferred 
through  points  four  and  five  of  Mesh  2  which  “donate”  values  to  points  IL  —  1  and  IL  of 
Mesh  1.  This  exchange  of  information  is  shown  with  arrows  at  each  point  in  the  schematic. 
Note  that  although  points  on  line  IL  —  2  (Mesh  1)  and  three  (Mesh  2)  are  also  collocated, 
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(a)  Swirl  velocity  at  T  =  12  along  horizontal 
centerline  with  several  schemes  coupled  with 
LOC  approach 


(b)  Vorticity  contours  at  three  instants  with 
compact  difference  scheme 


Figure  3:  Multidomain  Computation  with  Three-Point  Vertical  Overlap 


they  do  not  communicate  directly.  This  type  of  point  only  occurs  when  the  number  of 
overlap  points  is  odd  and  will  be  designated  the  noncommunicating  overlap  (NCO)  point: 
it  facilitates  the  detection  of  “drift”  between  the  two  solutions. 

Larger  and  smaller  overlap  regions  can  also  be  addressed  in  a  similar  manner  with  a 
minimum  of  at  least  two  points.  At  this  time  overlaps  with  only  an  odd  number  of  points 
have  been  investigated,  specifically  those  consisting  of  three,  five  and  seven  points.  The 
conclusions  derived  below  are  anticipated  to  carry  over  to  cases  with  even  number  of 
overlap  points. 

2.2.7  Three-point  overlap  results 

Figure  3(a)  exhibits  the  swirl  velocity  along  the  centerline  for  a  three-point  overlap  i.e.,  one 
where  each  domain  obtains  information  from  the  other  along  a  single  line.  Not  including 
the  overlap,  each  mesh  now  consists  of  30  x  30  points  with  the  NCO  occurring  at  X  =  6. 
multidomain  calculations  with  34443  and  66666  difference  schemes  are  plotted  with  theory 
and  the  single-domain  calculation,  which  is  the  best  achievable  result  for  this  purpose.  The 
differencing  schemes  in  the  multidomain  calculations  differ  substantially  in  formal  order  of 
accuracy,  the  first  34443  is  formally  fourth  order  accurate  by  the  results  of  Ref.  26  while 
the  second,  66666,  is  formally  sixth  order  accurate.  Thus,  the  large  swirl  velocity  errors  in 
the  multidomain  calculations  are  not  caused  by  the  difference  scheme.  Subsequent 
investigations  therefore  consider  only  the  34443  scheme.  Vorticity  contours  at  three 
time-instants,  shown  in  Figure  3(b)  manifest  a  large  degree  of  smearing  and  reduction  in 
peak  vorticity  values  as  the  disturbance  crosses  the  overlap.  Note  particularly  the  contours 
at  T  =  6  where  the  vortex  appears  to  elongate  as  it  passes  through  the  interface.  The  most 
likely  cause  of  error  is  the  lower  order  LOC  Liter  formulation  employed  in  each  domain 
near  the  interface. 
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(a)  Swirl  velocity  at  T  =  12  along  horizontal 
centerline 


0  10  20 


(b)  Vorticity  contours  at  three  instants  with 
compact  difference  scheme 


Figure  4:  Effect  of  Higher  Order  Formulas  With  three-Point  Vertical  Overlap 


To  test  this  hypothesis,  the  LOO  approach  is  replaced  with  higher  order  one  sided  filter 
formulas  developed  as  part  of  this  effort  and  presented  in  Ref.  9.  The  results  show 
significant  improvement  in  accuracy.  Figure  4(a)  displays  the  computed  swirl  velocity  along 
the  centerline  with  a  variety  of  higher  order  interface  treatments  together  with  the  baseline 
single-domain  calculation.  The  plot  indicates  that  raising  the  minimum  order  of  accuracy 
to  fourth  improves  performance  substantially  with  peak  error  in  swirl  velocity  of  only  6 
percent  compared  with  40  percent  for  the  second  order  LOG  method  (Figure  3(a)).  With  a 
minimum  sixth  order  filter,  FTO0'4  —  0  •  6  ■  6  •  6  •  8,  the  accuracy  is  comparable  to  the 
original  single-domain  calculation.  Further  increase  in  minimum  filter  order  of  accuracy  to 
eighth  and  tenth  order  has  no  impact  on  the  solution,  which  remains  stable.  Vorticity 
contours  with  the  FTO0'4  —  0  ■  8  •  8  •  8  •  8  scheme  are  shown  in  Figure  4(b)  and  compare  very 
well  with  those  obtained  with  the  single-domain  calculation  shown  in  Figure  1(a). 

2.2.8  Effect  of  Overlap  Size 

In  a  three-point  overlap,  each  domain  receives  data  at  an  end  point,  which  as  previously 
stated  is  not  filtered.  By  contrast,  with  larger  overlaps  (Figure  2),  each  mesh  receives  data 
at  several  near-boundary  points  which  are  set  to  values  obtained  from  adjacent  domains. 
The  filter  treatment  of  these  near-boundary  points  propagates  into  the  interior  due  to  the 
implicit  nature  of  the  Pade-type  formulas.  In  fact,  as  discussed  later,  the  LOG  approach 
can  actually  introduce  instabilities  at  domain  interfaces.  Based  on  heuristic  studies,  points 
receiving  data  from  adjacent  domains  should  preferably  be  left  unfiltered.  However,  higher 
order  one  sided  formulas  may  also  be  employed.  In  either  case,  the  minimum  order  of 
accuracy  increases  with  overlap  size,  and  consequently  interface  treatment  becomes  less 
critical.  This  is  demonstrated  again  by  plotting  swirl  velocity  along  the  vortex  centerline  in 
Figure  5  for  the  three,  five  and  seven  point  overlaps.  The  seven  point  overlap  computation, 
which  has  a  minimum  sixth  order  filter,  essentially  recovers  the  single  domain  calculation. 
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Figure  5:  Effect  of  Overlap  Size  on  Swirl  Velocity  Along  Centerline  at  T  =  12 


A  more  challenging  problem  consists  of  the  case  where  the  original  single-domain  mesh  is 
split  horizontally  so  that  the  top  part  of  the  vortex  is  always  in  Mesh  1  while  the  lower 
part  is  always  in  Mesh  2  (Figure  6(a)).  The  performance  of  several  schemes  with  a  five 
point  horizontal  overlap  is  shown  in  Figure  6(b).  The  minimum  fourth  order  (LOG) 
formulation  shows  considerably  higher  loss  of  accuracy  compared  to  the  five  point  vertical 
overlap  calculation  of  Figure  5.  Increasing  the  minimum  order  of  accuracy  brings  the 
results  successively  closer  to  the  single-domain  result.  There  is  little  difference  between  the 
minimum  eighth  and  tenth  order  formulations.  Vorticity  contours  obtained  with  the  higher 
order  scheme  exhibit  negligible  distortion  as  shown  in  Figure  6(a).  In  contrast,  results  with 
the  LOC  approach  (not  shown)  indicate  that  the  vortex  is  diffused  and  distorted  in  the 
direction  normal  to  the  interface. 


2.2.9  Curvilinear  Multidomain  Calculations 


To  demonstrate  the  use  of  the  interface  treatment  for  a  general  curvilinear  mesh,  consider 
the  wavy  grid  shown  in  Figure  7(a).  The  mesh  is  generated  with  an  extension  of  the 
formula  described  in  Ref.  23: 


xi,j 


Vi,j 


•Emin  H- 
Vmin  T  Ay0 


\i  -  1)  +  Axsm(nxW{j-^yo 
\j  -  1)  +  Aysm{  n^l)Axo 


+ 

+ 


i&x  \ 

IL—1 1 
jj’y  \ 
JL-1> 


(20) 
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(a)  Contour  plots  of  vorticity  with  FIO0'4  —  0  • 
0  •  8  •  8  •  8 


(b)  Swirl  velocity  at  T  =  12  along  interface 
boundary 


Figure  6:  multidomain  Computations  with  Five-Point  Horizontal  Overlap 


At  =  — —  /V?/  =  — — 
L^dj°  IL- l’^o  ,/L-l 

1  <  %  <  IL,  1  <  j  <  JL. 


(21) 


where  JL,  JL  denote  the  number  of  points  in  the  x  and  y  directions  and  Lx  =  xmax  —  xmin 
and  Ly  =  ymax  —  Vmin  are  the  corresponding  nominal  ranges.  The  mesh  of  Figure  7(a)  was 
generated  with  the  parameters  IL  =  60,  JL  =  30,  xmin  =  ymin  =  —6,  xmax  =  18,  y  max  6  ? 
Ax  =  1,  Ay  =  4,  nx  =  4  =  ny  =  4,  <px  =  <py  =  0.  The  resulting  mesh  has  a  maximum 
deviation  from  orthogonality  of  63.1  degrees  and  a  ratio  of  minimum  to  maximum  Jacobian 
of  2.1.  As  before,  the  34443  scheme,  is  retained  with  a  five  point  overlap,  v  components  of 
the  velocity  along  the  line  j  =  JL/ 2  at  T  =  12  computed  with  several  interface  filter 
formulations  are  shown  in  Figure  7(b)  together  with  the  single-domain  calculation.  Each 
scheme  shows  some  loss  of  accuracy  compared  with  the  equivalent  uniform  mesh  result  (see 
Figure  1(b)).  This  is  the  consequence  of  the  curvilinear  nature  of  the  mesh,  an  issue 
discussed  in  Ref.  21.  The  pertinent  point  here  however  is  that  even  on  this  wavy  mesh,  the 
higher  order  interface  treatments  are  stable  and  recover  the  single-domain  solution.  The 
discrepancy  in  peak  swirl  velocity  between  the  single-  and  the  two-domain  computation 
with  LTO0'4  —  0  ■  0  •  8  •  8  •  8-  is  less  than  1  percent.  Figure  7(c)  displays  contours  of  the 
perturbation  velocity  \J {u  —  u^)2  +  v 2  obtained  with  the  LTO0'4  —  0  •  0  •  8  •  8  •  8  filter.  The 
vortex  is  well-behaved  as  it  moves  through  the  interface.  The  contour  lines  are  continuous 
at  the  noncommunicating  overlap  point. 


2.2.10  Viscous  Multidomain  Calculations 

The  experience  gained  from  the  previous  inviscid  cases  is  now  applied  to  a  viscous  flow  A 
stringent  case  is  constructed  by  splitting  the  domain  horizontally,  in  such  a  manner  that 
the  interface  passes  through  the  boundary  layer  as  sketched  in  Figure  8(a).  The  implicit 
time-integration  algorithm  in  the  present  code  mandates  a  minimum  five-point  overlap. 
The  compact  fully  fourth  order  44444  difference  scheme  is  coupled  with  an  eighth  order 
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5-POINT  OVERLAP 


0  10  20 


(a)  Mesh 


(b)  a-component  of  velocity  at  T  =  12 


(c)  Perturbation  velocity  contours  at 
three  time-instants 

Figure  7:  Propagation  of  Vortex  on  Multidomain  Curvilinear  Mesh 

interior  filter  («/  =  0.4).  Since  the  latter  requires  a  9-point  stencil,  the  lower  domain  (#1) 
is  chosen  to  extend  a  total  of  9  points  normal  to  the  plate.  In  this  domain,  the  interior 
filter  formula  is  applied  only  at  point  five  since  all  others  points  require  near-boundary 
treatment.  The  noncommunicating  overlap  point,  7,  is  located  at  y/5  ~  0.57  at  the 
streamwise  section  corresponding  to  a  Reynolds  number  of  1.45  x  105.  The  outer  domain 
(^2)  extends  from  point  five  of  the  original  mesh  to  the  outer  boundary. 

Even  on  this  highly  stretched  mesh,  no  numerical  instabilities  are  detected  due  to  the 
introduction  of  the  interface.  The  success  of  the  higher  order  boundary  treatment  and  its 
superiority  over  the  LOC  approach  is  evident.  Close  examination  reveals  that  higher  order 
results  do  not  show  any  significant  “drift”  in  solutions  at  the  NCO.  Note  that  in  the  LOC 
results,  the  F80'3  —  0  •  2  •  4  •  6  •  8  filter  is  employed  instead  of  the  more  appropriate 
F80'3  —  0  ■  0  ■  4  •  6  •  8  suggested  earlier.  This  latter  choice  however  is  unstable  since  it  leaves 
point  two  unfiltered  due  to  present  code  limitations  requiring  a  “symmetric”  choice  of 
Liters. 

The  accuracy  and  robustness  of  the  boundary  filtering  formulation  is  further  examined  by 
considering  the  more  practical  case  of  vortex  shedding  past  a  circular  cylinder.  The  specific 
flow  conditions  are  freestream  Mach  number  M <*,  =  0.1  and  Reynolds  number  Rep  =  100, 
based  on  cylinder  diameter.  A  nonorthogonal  O-grid  of  size  155  x  201  is  employed  with  a 
nonuniform  spatial  distribution  in  both  the  radial  and  circumferential  directions.  All 
results  presented  are  obtained  with  the  implicit  Beam- Warming  time-marching  algorithm 
using  two  sub- iterations  and  a  time  step  IStU^/D  =  0.002. 

Since  this  flow  held  exhibits  spatial  periodicity  in  the  azimuthal  direction,  it  is  first 
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Implicit  time-integration 

Domain  1 

Boundary  layer 

Interface:  5  point  overlap 

l 

^^/^Domain  2 

V 

(a)  Schematic  (b)  Velocity  profiles  on  mesh  FPBL1 

Figure  8:  Multidomain  Flat-Plate  Boundary  Layer  Computation  with  Five-Point  Overlap 


computed  with  a  fully  periodic  formulation  of  the  compact  C4  scheme  and  the  eighth  order 
filter,  F 80'3.  Contours  of  the  instantaneous  vorticity  are  shown  in  Figure  9(a)  depicting  the 
Karman  vortex  street.  As  previously  shown  in  Ref.  23,  these  results  were  found  to  be  in 
excellent  agreement  with  available  highly-resolved  computations,27  in  terms  of  Strouhal 
number  and  lift / drag  coefficients.  In  order  to  evaluate  the  one  sided  high  order  filter 
formulation,  the  flow  was  recomputed  without  invoking  the  property  of  spatial  periodicity. 
Instead,  a  cut  was  introduced  in  the  azimuthal  direction  with  a  five  point  overlap,  as 
illustrated  in  the  schematic  of  Figure  9(b).  This  cut  was  deliberately  placed  in  the  middle 
of  the  vortex  street  to  provide  a  more  stringent  test  of  accuracy  and  stability.  Results 
computed  with  44444  and  a  minimum  sixth  order  filter  formulation  (F 80,3  —  0  ■  6  •  6  •  6  •  8) 
are  shown  in  Figure  9(c).  These  results  are  essentially  identical  to  the  fully  periodic 
solution  of  Figure  9(a).  By  contrast,  the  LOC  method  is  found  to  be  unstable.  Stability 
(and  accuracy)  can  be  recovered  by  heuristic  adjustments  of  near-boundary  values  of  af  in 
this  case,  it  is  necessary  to  decrease  this  coefficient  at  points  three  and  IL  —  2.  This 
application  again  demonstrates  the  superiority  and  elegance  of  the  one  sided,  higher  order 
approach  since  no  local  adjustment  of  the  filter  coefficient  is  required. 

2.3  Applications 

The  truly  multidisciplinary  capability  of  the  above  algorithm  is  highlighted  with  several 
examples  in  fluid  dynamics,  electromagnetics,  and  magnetogasdynamics.  Results  pertinent 
to  the  requirements  of  fluid-structure  interactions  may  be  found  in  Ref.  19. 

2.4  Fluid  Dynamics 

The  first  calculation  addresses  the  unsteady  simulation  of  spiral  vortex  breakdown  above  a 
slender  delta  wing  with  a  sweep  angle  of  75°  and  at  32°  angle  of  attack.  The  freestream 
Mach  number  is  0.2  and  the  Reynolds  number  based  on  centerline  chord  is  9200.  Analysis 
of  this  complex  phenomenon,  which  has  a  significant  impact  on  aircraft  control  and 
structural  integrity,28  remains  a  challenge,  partly  due  to  the  prohibitively  large 
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Arrows  show 
information  transfer 


(c)  Higher  order  one  sided  formulations  at 
overlap 


Figure  9:  Shedding  Past  a  Cylinder  at  Re  =  100 

computational  resources  required  for  accurate  studies.  An  extensive  computational  analysis 
with  the  E2  scheme,  including  a  mesh  resolution  study,  has  been  reported  in  Ref.  28  and  is 
employed  to  evaluate  the  accuracy  and  efficiency  gains  of  the  new  higher  order  scheme. 

The  delta  wing  is  embedded  in  an  H-H  topology,  highly  stretched  mesh  and  represents  a 
challenge  in  the  implementation  as  well  as  a  test  of  the  robustness  of  the  present  scheme. 
Figures  10  (a),  (b)  show  the  instantaneous  spanwise-component  of  vorticity  on  a  vertical 
plane  cutting  through  the  vortex  core  with  E2  for  two  different  grids.  Grid  2  contains 
double  the  resolution  of  Grid  1  in  the  core  region  and  was  shown  in  Ref.  28  to  be  adequate. 
At  this  angle  of  attack,  the  breakdown  location  is  highly  sensitive  to  details  of  the 
numerical  simulation.  Thus,  on  the  coarse  mesh,  the  E2  scheme  results  in  a  premature 
breakdown  of  the  vortex.  Figure  10  (c)  displays  results  with  the  C6  scheme  supplemented 
with  the  tenth  order  filter  on  the  same  mesh.  It  is  observed  that  the  breakdown  location  is 
similar  to  that  obtained  with  E2  on  the  fine  mesh.  Whereas  Figure  10  (a)  shows  relatively 
little  detail  with  only  two  smeared  concentrations  of  azimuthal  vorticity  in  the  vortex 
wake,  the  sixth-order  scheme  C6  with  a  tenth-order  filter,  Figure  10  (c),  exhibits  a  much 
richer  structure  corresponding  to  a  stronger,  tightly  wound  spiral  which  is  in  far  better 
agreement  with  the  grid-converged  pattern  of  Figure  10  (b).  An  estimate  based  on  these 
results  suggests  that  the  higher  order  method  yields  a  roughly  eight-fold  reduction  in  mesh 
points  and  sixteen-fold  reduction  in  computational  time. 

As  another  example,  the  method  has  been  employed  to  simulate  the  three-dimensional 
transition  of  a  forced,  finite  aspect  ratio,  plane  wall  jet.29  The  wall  jet  configuration 
considered  is  shown  in  Figure  11  (a).  The  main  parameters  governing  the  flow  are  the 
Reynolds  number,  the  disturbance  characteristics  at  the  jet  nozzle,  the  aspect  ratio  of  the 
channel,  and  the  length  of  the  wall.  In  the  present  study,  the  Reynolds  number,  based  on 
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Figure  10:  Application  to  Fluid  Dynamics:  Vorticity  Field  in  Vortex  Breakdown  over  Delta 

Wing 


Figure  11:  Application  to  Fluid  Dynamics:  Simulation  of  a  Forced  Transitional  Wall- Jet 
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exit  maximum  velocity  ( Umax )  and  nozzle  height  (h),  is  2150.  The  forcing  frequency  for 
this  case  is  200 Hz  while  its  amplitude  is  6  percent  of  the  maximum  centerline  velocity.  The 
aspect  ratio  of  the  channel  is  2 b/h  =  20.  The  mean  velocity  profile  in  the  normal  direction 
at  the  nozzle  exit  is  parabolic  and  corresponds  to  a  fully-developed  laminar  channel  flow  in 
that  direction. 

An  example  of  the  ability  of  the  present  high-order  scheme  to  discern  the  fine-scale 
breakdown  to  turbulence  is  shown  in  Figure  11  (b)  which  depicts  a  3-D  representation  of 
the  structure  of  shear-layer  rollers  in  terms  of  an  iso-surface  of  vorticity  magnitude.  The 
transition  process  begins  with  the  formation  of  shear-layer  and  wall  vortex  pairs  which,  due 
to  the  energetic  forcing,  appear  close  to  the  nozzle  exit  and  are  phased-locked  for  a  short 
distance  downstream.  In  the  course  of  their  spanwise  evolution,  the  rollers  are  first  split 
into  a  double-helical  structure,  which  is  clearly  discernable  near  the  sidewalls.  This  feature 
propagates  toward  the  center  while  also  expanding  in  the  radial  direction.  The  spiral 
branches  are  wound  in  a  sense  opposite  to  that  of  the  swirl  direction  of  the  vortex, 
consistent  with  the  direction  of  the  induced  axial  flow  which  exists  within  the  vortex  core. 
This  vortex  branching  and  helical  twisting  spreads  rapidly  through  self-induction  effects, 
and  eventually  reaches  the  symmetry  plane  where  the  vorticity  magnitude  within  the 
vortex  core  is  drastically  diminished  (hence  the  apparent  break  in  the  iso-surface,  region  ‘y’ 
in  Figure  11  (b).  Additional  details  of  the  simulation  and  the  physics,  including  impact  of 
forcing  parameters  and  endwall  conditions,  may  be  found  in  Ref.  29. 

2.5  Wave  Propagation  —  Electromagnetics 

The  low  dispersion  error  characteristic  of  compact-difference  schemes  is  an  attractive 
property  in  the  simulation  of  wave  propagation  phenomena  associated  with  acoustic  and 
electromagnetic  scattering.  In  order  to  demonstrate  the  capability  of  the  present  numerical 
approach  to  treat  such  phenomena,  consider  the  problem  of  electromagnetic  scattering  by  a 
sphere  of  diameter  a  =  1  as  an  illustrative  example  demonstrating  the  capability  of  the 
current  scheme  to  solve  Maxwell’s  equations.  The  domain  of  computation  is  truncated  at 
an  outer  spherical  boundary  placed  at  r  =  2.5  and  is  discretized  by  a  spherical  49  x  49  x  99 
mesh  with  uniform  variation  of  dr,  dd  and  d<p.  The  scattered  field  formulation  is  employed 
with  a  transverse-magnetic  (TM)  Gaussian  pulse  for  illumination,  and  the  radar 
cross-section  (RCS)  is  obtained  by  sampling  the  response  at  various  frequencies.  The 
bistatic  RCS  profiles  at  ak  =  4  and  ak  =  7,  where  k  is  the  wave  number,  are  shown  in 
Figure  12  for  the  0  =  0  longitudinal  plane.  At  these  two  values  of  ak,  the  number  of  points 
per  wave  (PPW)  in  the  radial  direction  varies  from  19  PPW  at  ak  =  4  to  11  PPW  at 
ak  =  7  while  the  outer  boundary  placement  varies  from  2.55A  to  4.46A  respectively.  The 
results  are  seen  to  be  in  very  good  agreement  with  the  Alie  theory  and  demonstrate 
considerable  advantage  over  traditional  methods  which  require  over  25  PPW  for  adequate 
resolution.30 
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Figure  12:  Application  to  Electromagnetics:  Radar  Cross  Section  of  a  Sphere  at  ak  =  4.0 

and  ak  =  7.0 


2.6  Magnetogasdynamics 

There  is  presently  considerable  interest  in  exploring  the  possibility  of  utilizing 
electromagnetic  forces  to  control  conducting  fluid  flow,  particularly  in  the  hypersonic 
regime,  to  alleviate  the  stringent  environments  encountered  and  to  enhance  the  efficiency  of 
propulsion.  From  a  numerical  standpoint,  the  problem  requires  careful  simulation  of  the 
interaction  between  electromagnetic,  pressure  and  viscous  forces  described  by  the 
magnetogasdynamic  equations.  The  ability  of  the  present  scheme  to  solve  these  equations 
has  been  explored  extensively  in  Refs.  31,  32  and  17.  As  an  example,  consider  the  flow 
through  a  channel,  depicted  in  Figure  13  (a)  together  with  a  simplified  representation  of 
the  circuit  completing  the  current  path.  The  pertinent  parameter  in  such  flows  is  the 
Hartmann  number,  Ha  =  \/ ReReaRb  =  BrefL^J crref  / /aref  which  is  a  measure  of  the  ratio  of 
magnetic  to  viscous  forces.  A  uniform  magnetic  held  is  imposed  in  the  ^/-direction  (By  =  1) 
and  an  electric  held,  measured  by  the  nondimensional  parameter  K  =  —Ez/(UrefBref),  is 
imposed  across  the  channel.  If  end  effects  are  ignored,  the  fully  developed  velocity  profile 
can  be  obtained  by  analytic  means33  and  can  then  be  employed  for  validation  purposes. 

The  problem  is  modeled  by  considering  a  finite  streamwise  length  of  the  channel  as  shown 
in  Figure  13  (b).  The  traditional  zero  normal  pressure  gradient  boundary  conditions  are 
modified  to  include  the  term  due  to  magnetic  pressure  as  described  in  Ref.  32. 

Figure  13  (c)  exhibits  the  velocity  prohle  obtained  with  C4  and  the  eighth-order  filter  for 
K  =  2  and  Ha  =  100  (Re  =  100,  Rb  =  10,  Rea  =  10).  The  velocity  prohle  without  the  B 
held  is  also  shown  together  with  the  theoretical  result,  the  equation  for  which  may  be 
found  in  Ref.  33.  Three  different  meshes  are  considered  consisting  of  30  x  50,  50  x  80  and 
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(a)  Overall  flow  configuration  (b)  Problem  setup  (c)  Velocity  profiles 

Figure  13:  Application  to  Magnetogasdynamics:  Control  of  Flow  in  a  Channel  with 

Transverse  Magnetic  Field 


50  x  110  points.  The  finest  mesh  contains  about  35  points  in  the  Hartmann  layers  near 
walls,  ft  is  evident  that  the  two  coarser  meshes  exhibit  error  near  the  edge  of  the 
Hartmann  layer  (see  insert  of  Figure  13  (c)).  In  comparison  with  the  case  where  no  B  held 
is  applied,  the  effect  of  the  magnetic  held  is  to  hatten  the  velocity  as  expected  since 
Lorentz  forces  diminish  shear.  Only  near  the  channel  walls,  where  velocity  gradients  are 
high,  are  viscous  stresses  of  the  same  order  as  the  Maxwell  stresses.  The  computational 
efficiency  advantage  from  the  present  higher  order  method  for  this  case  is  similar  to  that 
reported  earlier  in  the  context  of  vortex  breakdown  over  the  delta  wing. 


21 


3  High-temperature  molecular  kinetics  methodology 


The  presence  of  shock  waves  in  high  speed  flow  of  an  air  mixture  presents  considerable 
difficulties  for  accurate  numerical  simulation  of  the  flowfield.  The  shock  wave  redistributes 
the  high  kinetic  energy  of  the  oncoming  flow  into  various  internal  energy  modes,  which 
relax  relatively  slowly,  leading  to  significant  chemical  and  thermal  nonequilibrium  in  the 
stagnation  region.  In  the  gas  kinetic  description,  intermolecular  collisions  change  the 
translational,  rotational,  vibrational,  and  electronic  energies  of  the  collision  partners.  The 
probabilities  or  effective  cross  sections  of  these  elementary  processes  differ  significantly, 
giving  rise  to  widely  separate  relaxation  times  for  the  internal  modes.  Thus  it  becomes 
important  to  account  for  the  rates  of  relaxation  processes  to  predict  the  nonequilibrium 
behavior  of  these  kinds  of  flows.  Equilibrium  is  established  relatively  fast  for  the 
translational  degree  of  freedom  compared  to  internal  degrees  of  freedom.  The  equilibrium 
Maxwellian  distribution  is  established  as  a  result  of  the  exchange  of  momentum  and  kinetic 
energy  among  particles.  The  relaxation  time  for  establishing  a  Maxwellian  distribution  in 
components  of  air  is  of  the  order  of  the  average  time  between  gaskinetic  collisions. 

7~trans  ~  7" gas 

Here  /  is  the  gaskinetic  mean  free  path,  v  is  the  average  particle  velocity,  N  is  the  particle 
number  density,  and  crgas  is  the  gaskinetic  collision  cross  section.  In  air  at  standard 
conditions,  l  ~  6  x  ICE6  cm  and  Ttrans  ~  10“10  sec.  Usually  the  gaskinetic  times  are  small 
in  comparison  with  the  flow  times  over  which  appreciable  changes  in  the  macroscopic 
parameters  of  the  gas,  density  or  energy,  take  place.  When  these  conditions  are  satisfied  it 
is  possible  to  assign  at  every  instant  of  time  a  “translational”  temperature,  which 
characterizes  the  average  kinetic  energy  of  translational  motion  of  the  particles.  However, 
when  the  velocity  distribution  function  of  the  particles  is  anisotropic,  as  is  the  case  in  low 
density,  hypersonic  nonequilibrium  flows,  multiple  translational  temperatures,  one  in  the 
direction  of  flow  and  the  other  in  the  direction  perpendicular  to  the  flow,  need  to  be 
delineated. 

For  hypersonic  flow  conditions,  besides  the  nonequilibrium  translation,  the  internal  energy 
modes  get  excited  and  it  becomes  important  to  model  the  rotational  and  vibrational 
energy  relaxation  and  the  coupling  of  those  modes  to  the  dissociation  process.  The 
influence  of  internal  energy  excitation  on  dissociation  is  an  area  of  active  research  and  is 
addressed  in  the  present  report.  Early  research  in  dissociation  kinetics  revealed  two 
divergent  theoretical  views  on  the  treatment  of  dissociation.  Qualitatively  stated,  one  is 
the  weak  bias  mechanism34  where  dissociation  proceeds  from  the  entire  vibrational 
manifold,  the  dissociation  from  the  low  levels  driven  by  rotation.  The  other  is  the  strong 
bias  mechanism35  based  on  the  assumption  that  dissociation  requires  vibrational  excitation 
to  within  kT  of  the  dissociation  limit.  Dissociation  is  a  reactive  process  with  an  activation 
energy  typically  over  an  order  of  magnitude  greater  than  the  energy  associated  with  V-T 
exchange.  Nonreactive  kinetic  processes  tend  to  bring  about  an  equilibrium  distribution  in 
the  vibrational  manifold,  whereas,  the  dissociation  process  perturbs  it.36  In  the  upper 
quantum  energy  levels,  the  vibrational-translational  (V-T)  kinetic  rate  processes  dominate 
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Figure  14:  Schematic  Showing  Multiquantum  Transitions  in  Vibrational  States  and 

Dissociation 


the  energy  transfer.  In  some  circumstances,  the  vibrational-vibrational  (V-V)  kinetic  rate 
processes  can  play  an  important  role  by  accelerating  the  reaction  rates.3'  The  other  issue 
in  the  afore-mentioned  vibration-dissociation  coupling  is  that  of  population  depletion38  or 
enhancement39  in  the  vibrational  manifold  as  a  result  of  the  perturbation  caused  by  the 
dissociation  process.  A  schematic  of  the  vibrational  potential  well  depicting  the  various 
processes  is  shown  in  Figure  14. 

When  dissociation  takes  place,  the  losses  from  the  vibrational  manifold  reduce  the 
vibrational  populations  below  the  equilibrium  values,  resulting  in  a  decrease  in  the 
dissociation  rate.  In  the  work  of  Josyula  and  Bailey,40  the  generalized  depletion  equations, 
considering  the  state-to-state  kinetics  of  dissociating  nitrogen,  were  solved  to  predict  the 
extent  of  population  depletion  in  the  vibrational  manifold.  The  relative  importance  of  V-V 
exchanges  was  quantified  for  dissociation  from  last  quantum  level.  However,  dissociation 
can  proceed  from  any  quantum  level.  Recent  calculations  based  on  the  Forced  Harmonic 
Oscillator  (FHO)  model  suggest  that  for  a  small  portion  of  the  shock  layer  behind  the 
blunt  body,  where  the  vibrational  temperature  is  less  than  half  the  translational 
temperature,  the  dissociation  proceeds  from  the  low  vibrational  levels.41  For  the  region  of 
the  shock  wave  where  the  vibrational  temperature  is  nearly  equal  to  the  translational 
temperature,  there  is  only  a  weak  vibrational  bias  and  dissociation  proceeds  from  all  the 
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levels.  However,  these  state-specific  predictions  have  yet  to  be  verified  in  an  experiment. 

The  V-T  and  V-V  transfers  can  occur  by  multiquantum  energy  exchanges  in  the 
vibrational  manifold.  Multiquantum  vibrational  energy  exchanges  in  highly  excited 
molecules  for  the  CO  jC02  gas  dynamic  laser  were  simulated  by  several  researchers.42-45 
However,  there  has  been  little  work  on  multiquantum  energy  exchanges  for  hypersonic  air 
flow  past  blunt  bodies.  Multiquantum  V-T  and  V-V  transition  modeling  performed  in 
Ref.  46  using  Forced  Harmonic  Oscillator  (FHO)  model47,48  concluded  that  dissociation  of 
nitrogen  in  stagnation  regions  occurs  by  multiquantum  transitions  from  the  vibrational 
levels  below  u=10  and  not  via  ladder  climbing  process. 

The  present  study  simulates  nonequilibrium  relaxation  in  high  speed  flows  considering 
multitranslational  and  rotational  temperature  modeling.  Vibrational  relaxation  and 
vibrational-dissociation  coupling  are  examined  in  the  context  of  a  master  equation 
approach  under  different  bias  conditions. 


3.1  Theoretical  Model 


The  Boltzmann  equation  expresses  the  behavior  of  many-particle  kinetic  system  in  terms  of 
the  evolution  equation  for  the  single  particle  gas  distribution  function.  The  simplified 
collision  model  given  by  the  Bhatnagar  Gross  Krook  (BGK)  model  is  formulated  as 


d£ 

dt 
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d  f  feq  ~  f 
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where  /,  the  distribution  function,  gives  the  number  density  of  molecules  at  position  x  and 
velocity  u  at  time  t.  The  left  hand  side  of  the  above  equation  represents  the  free  streaming 
of  molecules  in  space,  and  the  right  side  denotes  the  collision  term.  If  the  distribution 
function  /  is  known,  macroscopic  variables  of  the  mass,  momentum,  energy  and  stress  can 
be  obtained  by  integration  over  the  moments  of  molecular  velocity.  In  the  BGK  model,  the 
collision  operator  involves  a  simple  relaxation  to  a  local  equilibrium  distribution  function 
feq  with  a  characteristic  time  scale  r.  The  BGK  model  was  proposed  to  describe  the 
essential  physics  of  molecular  interactions,  with  r  chosen  as  the  molecular  collision  time. 
Based  on  the  above  BGK  model,  the  Navier-Stokes  equations  can  be  obtained  with  the 
Chapman- Enskog  expansion  truncated  to  the  lst-order, 


/  =  n  +  Knf\ 

=  feq -r(dfeq/dt  +  u-dfeq/dx).  (24) 

where,  feq  denotes  the  contribution  of  the  equilibrium  part  and  Kn  f\  the  nonequilibrium 
component.  With  the  modification  of  the  particle  collision  time,  the  validity  of  the  kinetic 
model  is  extended  beyond  that  of  the  Navier-Stokes  equations.  An  earlier  attempt  based 
on  this  model  was  successful  for  shock  structure  calculations  of  the  monatomic  gas  with  a 
single  translational  temperature.49 

With  the  above  kinetic  model,  the  Navier-Stokes  equations  in  mass-averaged  velocity  form 
can  be  derived, 

^(p)  + V- (pu)  =0  (25) 
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Further,  to  include  the  internal  energy  contribution  dne  to  the  inelastic  collisions  in  the 
vibrational  energy  mode,  the  above  mass  conservation  equation  (Eq.  25)  can  be  written  to 
include  the  source  term,  ion, 


d_ 

dt 


(. Pv )  +  V  ■  (p„u)  =  u, 


v=0,l,... 


(28) 


In  the  continuum  regime,  the  Navier-Stokes  equations  could  then  be  used  to  simulate  a 
flow  in  thermal  disequilibrium.  Equations  25  to  27  describe  the  conservation  of  mass, 
momentum  and  energy  in  the  flowfields  of  interest.  Eqn.  26  and  27  represent  the 
conservation  of  total  momentum  and  energy,  respectively.  A  microscopic  kinetic  approach 
was  taken  (Eqn.  28)  by  treating  the  molecule  as  anharmonic  oscillator,  calculating  the 
state  populations  using  the  master  equation.  In  the  treatment  of  vibrational  energy  for  the 
diatomic  species  in  the  master  equation  code,  a  separate  vibrational  conservation  equation 
is  not  necessary  as  the  vibrational  energy  can  be  calculated  for  each  quantum  level. 

The  above  equations  can  be  derived  from  the  kinetic  theory,  if  the  internal  energy  states  at 
different  quantum  levels  are  considered  as  different  species.  Modern  computational  fluid 
dynamics  method  consists  of  solving  the  continuum  equations  for  updating  the  total  mass, 
momentum,  and  energy  first  through  the  calculation  of  the  numerical  fluxes  at  each  cell 
interface,  the  source  term  given  by  the  vibrational  energy  master  equations,  J^( pv )  =  c uv 
can  then  be  updated  inside  each  cell.  The  advantage  of  the  above  continuum  formulation 
(Equations  25  to  27  and  Eq.  28)  is  that  nonequilibrium  internal  energy  effects,  e.g.  those 
due  to  the  vibrational  energy  mode,  can  be  described  in  a  detailed  manner. 

However,  the  nonequilibrium  effects  encountered  in  the  transition  to  the  rarefied  flows, 
such  as  translational  and  rotational  nonequilibrium,  cannot  be  satisfactorily  described  in 
the  continuum  model,  hence,  a  kinetic  description,  Eq.  23,  becomes  attractive.  In  the 
following  section,  the  modification  to  Eq.  23  to  formulate  a  multiple  translational 
temperature  model  is  presented. 


3.2  Multitranslational  Temperature  Modeling 

For  a  gas  with  different  translational  temperatures  in  the  x,  y,  and  z  directions,  an 
anisotropic  distribution  is  assumed  with  equilibrium  state  prevailing  in  a  given  direction, 
and  for  a  one  dimensional  flow,  written  as: 

9  =  p(  — )1/2(— )1/2(— )1/2exp[-A:r(n  -  U)2  -  Xyv2  -  A zw2},  (29) 

7T  7T  7T 

where  p  is  the  density,  U  is  the  macroscopic  velocity  in  the  x-direction,  and  (u,  v,  w)  the 
components  of  the  particle  velocity  in  the  x,  y  and  z  directions,  respectively.  The 
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parameter  A  is  related  to  the  gas  temperature,  given  by  Xx  =  M./2kTx,  \y  =  M./2kTy  and 
=  M./2kTz.  For  the  multitranslational  temperature  formulation  of  a  one  dimensional 
flow  considered  in  the  present  study,  we  can  assume  that  Ty  —  Tz,  such  that,  Xy  =  Xz.  The 
establishment  of  the  above  ellipsoid  Maxwell  distribution  is  assumed  to  be  a  consequence  of 
particle  collisions.  Over  a  long  period  of  time,  the  above  state  g  will  further  approach  an 
overall  equilibrium  state  given  by  g , 

Xeq 

g  =  p{— )3/2  exp[— Ae9((w  -  U )2  +  u2  +  w2)\.  (30) 

71 

If  the  mass,  momentum  and  energy  is  conserved  during  the  particle  collision  process,  the 
relationship  between  an  averaged  equilibrium  temperature  Teq  (  Xeq  =  M./2kTeq  )  and  the 
individual  temperatures  Tx,Ty,  and  Tz  is  given  by, 

Teq  —  \{TX  +  2Ty),  (31) 

where  the  assumption  Ty  =  Tz  for  the  one  dimensional  flow.  The  above  process  from  g  to  g 
can  be  described  by  the  gas-kinetic  model 

gt  +  ugx  =  (g  -  g)/r,  (32) 

which  is  like  the  BGK  model.50  From  the  above  kinetic  equation,  besides  the  Navier-Stokes 
equations  for  the  mass,  momentum  and  total  energy,  an  additional  thermal  energy 
equations  for  the  particle  random  motion  in  the  y  and  z  directions  is  obtained.  For  the  ID 
case,  the  thermal  energy  equation  becomes 

(pTy),  +  (pUTy)x  =  ^(Tx  -  T„),  (33) 

where  R  is  the  gas  constant  and  r  is  the  particle  collision  time.  The  two  temperature, 
continuum  Navier-Stokes  formulation  given  by  Eq.  32  resembles  the  hydrodynamical 
multiple  translational  temperature  equations  of  Candler  et.  al.51  The  solutions  of  the  shock 
structure  by  solving  Eq.  32  using  the  gas-kinetic  scheme52  has  similar  behavior  as  the 
results  of  multitemperature  Navier-Stokes  results  in  Ref.  51.  The  inclusion  of  two 
translational  temperature  as  in  the  above  model,  results  in  the  shock  structure  solutions 
similar  to  those  of  Ref.  51  which  are  thinner  with  a  nonphysical,  sharp  jump  in  the 
calculated  density  profile  inside  the  shock  layer,  inconsistent  with  the  DSMC  solutions. 

The  discrepancy  is  more  pronounced  for  the  high  Mach  number  cases.  Even  though  the 
equation  was  derived  from  the  Boltzmann  equation,  the  model  given  by  Eq.  32  to  derive 
the  macroscopic  multiple  temperature  hydrodynamical  model  in  Ref.  51  is  intrinsically 
flawed.  Unfortunately,  most  multitemperature  kinetic  models  in  the  literature  are  based  on 
the  same  assumption  of  an  ellipsoid  Maxwellian  as  in  Ref.  51. 

In  order  to  improve  the  previously  mentioned  kinetic  model  for  the  translational 
nonequilibrium  to  overcome  the  erroneous  shock  solutions,  we  reconstructed  the  particle 
collision  process.  Starting  from  the  gas  distribution  function  /  for  a  monatomic  gas,  the 
particle  collisions  drive  f  to  g  and  g:  instead  of  from  g  to  g  alone.  During  the  course  of 
particle  collisions,  it  is  difficult  to  distinguish  the  process  from  f  to  g  or  from  /  to  g.  In 
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terms  of  the  particle  collision  time  for  translational  nonequilibrium,  they  can  be  assumed 
to  be  the  same.  Note,  this  should  not  be  confused  with  the  nonequilibrium  internal  energy 
distribution  in  a  diatomic  gas,  where  the  collision  time  for  the  rotational  relaxation  may  be 
much  larger  than  the  translational  nonequilibrium  relaxation.  Therefore,  for  the 
monatomic  gas  we  construct  a  generalized  BGK  model  for  the  translational  nonequilibrium, 

ft  +  «/.  =  h—  +  —},  (34) 

2  T  T 

where  /  is  the  real  gas  distribution  function,  and  g  and  g  are  the  corresponding 
Maxwellians.  Note  that  the  particle  collision  time  r  is  related  to  the  local  dynamical 
viscosity  p  and  pressure  p ,  i.e.,  r  =  p/p.  The  relation  between  mass  p,  momentum  pU, 
total  energy  pE ,  and  the  thermal  energy  pEy_z  in  the  y  and  2  directions  are 

P 

pU 
pE 

pEy—Z 

where  xf  has  the  components 

=  (1  ,u,^(u2 +  v2 +w2),^(v2  +w2))T.  (36) 

Since  only  mass,  momentum  and  total  energy  are  conserved  during  particle  collisions,  / 
and  g  satisfy  the  condition 

J ~f)  +  (9  ~  fMadudvdw  =  S  =  (0,0,0,  s)T,  a  =  1,2,  3, 4.  (37) 

The  source  term  s  can  be  still  be  modeled  as  s  =  pR(Tx  —  Ty)/3r. 

For  any  kinetic  model  for  multiple  translational  temperature  nonequilibrium,  one  basic 
requirement  is  that  the  averaged  temperature  Teq  =  \{TX  +  2 Ty)  in  the  ID  case  should  be 
the  same  temperature  as  the  one  temperature  gas  kinetic  model.  If  we  take  moments 
(u  —  U)2,v2  and  w 2  on  the  model  given  by  Eq.  34,  we  can  obtain  the  temperature  evolution 
equations  for  individual  direction.  When  we  add  them  together  and  use  the  condition 
3 T  =  Tx  +  Ty  +  Tz,  a  single  equation  for  the  averaged  temperature  T  can  be  obtained, 

(pT),  +  ( pUT)x  =  (pTeq  -  pT)/T, 

which  is  the  same  temperature  evolution  equation  from  a  single  temperature  BGK  model, 
i.e.  ft,  +  ufx  =  (g  —  f)/r.  So,  the  averaged  temperature  from  the  two-temperature  kinetic 
model,  Eq.  34,  is  the  same  as  that  from  a  single  temperature  model.  The  model  given  by 
Eq.  34  is,  thus,  shown  to  be  a  valid  extension  to  recover  the  multiple  temperature 
nonequilibrium  from  the  original  BGK  model. 

In  order  to  test  the  above  model  (Eq.  34),  the  numerical  method  developed  in  Ref.  52  is 
used  to  solve  Eq.  34.  For  a  finite  volume  method, 

1  /'At 

Wp1  =  w;  +  —  |  (F^1/2(t)-Fj+1/2(t))dt  +  Spt,  (38) 
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where  WJ1  is  the  cell  averaged  mass,  momentum,  total  energy,  and  the  thermal  energy  in  y 
and  z  directions,  and  Fj+ 1/2  is  the  corresponding  fluxes  at  a  cell  interface  by  solving  Eq.  34. 
Note  that  At  is  the  time  step  At  =  tn+1  —  tn ,  and  S™  is  the  source  term  for  the  thermal 
energy.  The  evaluation  of  the  fluxes  is  based  on  the  gas  distribution  function  /  at  a  cell 
interface.  For  the  shock  structure  calculation,  due  to  the  smoothness  inside  the  shock  layer, 
the  gas  distribution  function  at  a  cell  interface  is  constructed  using  the  Chapman-Enskog 
expansion  and  has  the  form 


11  1 

/  =  2  (j  +  g)  -  2Tlft  +  ft  +  '“(ft  +  ft)]  +  5*  (ft  +  ft) 


(39) 


where  the  gx  and  gx  can  be  obtained  from  the  gradients  of  macroscopic  variables,  and  gt 
and  gt  can  be  evaluated  from  the  compatibility  conditions 

J  (gt  +  ugx )  -0Q  du  dv  dw  =  0  and  J  (gt  +  ugx)dudvdw  =  0. 

For  more  details,  see  Ref.  52.  Even  with  the  ability  to  recover  two  translational 
temperatures,  the  nonequilibrium  distribution  (Eq.  39)  is  truncated  to  the  Navier-Stokes 
order  and  it  will  not  be  accurate  in  the  description  of  rarefied  flow.  The  anomalous  density 
jump  given  by  the  model  of  Eq.  32(found  in  Ref.  51)  does  not  appear  in  the  solutions  given 
by  the  present  formulation.  However,  the  shock  thickness  is  still  thinner  in  comparison 
with  experimental  measurements  and  the  DSMC  solutions.  Theoretically,  it  seems  that  we 
should  continue  this  expansion  in  Eq.  39  to  go  to  Burnett  and  Super-Burnett  orders. 
Unfortunately,  the  success  from  the  higher  order  expansion  is  limited,  and  there  is  no 
evidence  that  the  successive  expansions  converge.  Additionally,  we  can  replace  the  collision 
time  r  in  Eq.(39)  with  a  generalized  one,  t*,  and  assume  that  there  is  a  closed  form 
solution  of  the  BGK  model.  We  can  obtain  a  relation  between  r*  and  r.  In  the  current 
study,  the  new  collision  model  (Eq.  34)  is  used  with  BGK-Xu  method49  to  establish  a 
generalized  particle  collision  time  r*, 


r*  =  r/(l  +  r  <  D2g  >  /  <  Dg  >) 


(40) 


where  D  =  dt  +  udx  and  <  ...  >=  —  U)2dudvdw.  Here  r  depends  on  the 

macroscopic  variables  through  the  relation  rp  =  /j,00(T/T00)U)  ,  and  the  generalized  t*  will 
depend  not  only  on  the  macroscopic  variables,  but  also  on  their  gradients.  Therefore,  the 
evaluation  of  r*  represents  a  set  of  generalized  constitutive  relations  for  the  calculation  of 
the  translational  nonequilibrium  flow.  In  order  to  simulate  the  flow  with  a  realistic  Prandtl 
number,  a  modification  of  the  heat  flux  in  the  energy  transport  at  a  cell  interface,  as  used 
in  Ref.  52,  is  implemented  in  the  current  study. 


3.3  Rotational  Temperature  Modeling 

The  model  given  by  Eq.  34  can  be  extended  to  the  rotational  nonequilibrium  as  well.  As 
shown  above,  the  average  temperature  of  the  multiple  translational  model  is  the  same  as 
the  temperature  for  one  translational  temperature  model.  Therefore,  we  will  use  a  single 
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Figure  15:  Multiquantum  V-T  Transition  Rates  for  five  Quantum  Levels  at  1000  K  for  the 
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temperature  for  the  translational  degree  of  freedom  in  what  follows.  With  the  inclusion  of 
rotational  degree  of  freedom,  the  equilibrium  state  can  be  written  as 

9  =  9(Lp(L)  exp[-A,((«  -  U)2  +  v2  +  w2)  -  Ar£2].  (41) 
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Due  to  the  particle  collisions,  the  above  ellipsoid  distribution  will  approach  equilibrium, 

\eq 

g  =  p{— )5/2  exp[-A£9((u  -  U )2  +  v2  +  w2  +  f2)].  (42) 

7T 

Similar  to  the  model  given  by  Eq.  34,  we  can  construct  a  kinetic  model  for  translational 
and  rotational  nonequilibrium, 

/.  +  <*/„  =  h—  +  — ].  (43) 

2  Tt  Tr 

where,  clue  to  the  different  time  scales  of  translational  and  rotational  modes  to  achieve  the 
equilibrium  state,  two  different  collision  times  rt  and  rr  are  used  in  the  present  study  for 
the  condition  rr  »  Tt.  The  relationship  of  rr  =  ZrTt  was  invoked,  with  the  value  of  Zr 
ranging53  between  three  and  five  for  the  flows  considered  in  the  present  study. 

The  distribution  function  /  can  approach  g  rapidly,  then  slowly  from  g  to  g.  Thus,  the 
above  model  Eq.  43  can  be  further  simplified  as  shown  below  to  calculate  the  numerical 
fluxes, 

ft  +  ufx  =  (g  -  f)/rt,  (44) 
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where  the  fast  transition  part  is  captured  in  the  above  model.  For  the  rotational 
nonequilibrium  flow,  a  finite  volume  scheme  similar  to  Eq.  38  was  used  with  the  numerical 
fluxes  based  on  the  above  simplified  model  Eq.  44. 


3.4  Multiquantum  Energy  Exchanges  in  the  Vibrational 
Manifold 


The  conservation  Eqn.  28  is  written  for  the  mass  density  in  quantum  level  v.  The  source 
term  ujv  derived  from  the  vibrational  master  equations  includes  the  relevant  energy 
exchange  processes  consisting  of  the  V-T,  V-V  and  dissociation.  The  mass  density  of  the 
molecular  species  is  the  sum  of  the  corresponding  state  densities  in  the  vibrational  levels. 

vlast 

pv  (45) 

v=0 

The  symbolic  equations  governing  the  V-T  transitions  responsible  for  the  variation  of  the 
particles  distributed  in  the  Vth  vibrational  level  are: 

N2(v)  +  No  n  N2(v')  +  N2  (46) 

the  equations  governing  the  V-V  process  are: 

N2(v)  +  N2(w)  fl  N2(v')  +  N2(w')  (47) 

and  the  equations  governing  the  dissociation  process  for  the  Ladder  model  shown  here  for 
reference  (the  following  section  presents  the  weak  bias  model  as  well): 

N 2  (viast)  +  N2  — >  2 N  +  N2  (48) 


The  kinetics  of  the  particle  exchanges  among  the  quantum  states  are  simulated  using  the 
vibrational  master  equation,  the  population  distributions  are  calculated  with:54 


iov  =  -j- {Y^vtW  ->  v)pc’p  -  kVT(v  -►  v')pvp]  + 

J  v' 

y  [rvviv^w1  — >  v,w)pv>pw‘  -  rvv(v,w  — >  v\w’)pvpw\} 


(49) 


The  V-T  process  is  associated  with  the  rate  coefficient  kyr  where  the  molecule  loses  or 
gains  a  vibrational  quantum.  The  de-excitation  rate  from  v'  to  v  is  denoted  by 
kvrW  v),  the  inverse  collision  from  v  — >  v'  by  kyriy  — >  v').  The  multiquantum 
transition  rates  are  from  the  work  of  Adamovich,  et  al.48  shown  in  Figure  15,  and  16,  are 
for  1,000  K  (5  quantum)  and  10,000  K  (5  quantum).  At  1,000  K,  the  single  quantum  V-T 
transition  rates  of  Capitelli,  et  al.  match,  whereas,  at  10,000  K  the  Capitelli,  et  al.  rates 
are  higher  than  the  FffO  rates,  with  the  maximum  difference  an  order  of  magnitude  for  the 
last  level.  On  the  other  hand,  when  considering  V-V  exchanges,  the  initial  and  final 
vibrational  states  of  each  collision  partner  must  be  identified;  thus  the  transition  rate  from 
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Multiquantum  (5)  rates  at  10000K 


Figure  16:  Multiquantum  V-T  Transition  Rates  for  Five  Quantum  Levels  at  10000K  for 

the  1V2  Molecule 


v'  to  v  and  w'  to  w  is  given  by  ryV(y' ,w'  — >  v,w).  Consistency  of  the  rate  coefficient  with 
the  principle  of  detailed  balance  is  enforced.  Single  quantum  V-V  transfer  rates 
implemented  in  the  present  study  were  taken  from  Doroshenko,  et  al.55 

The  vibrational  energy  of  the  diatomic  molecules,  treated  as  anharmonic  oscillators,  is 
given  by, 

v°  o 

eVib  =  £  ~ev  (50) 

where  the  index  v  enumerates  the  vibrational  quantum  level.  In  this  equation,  pv/ p  is  the 
fractional  population  of  the  Vth  vibrational  level  and  ev  the  quantum  level  energy  given  by 
the  third-order  approximating  formula: 

11  1 

jj~c  =  ue(v  +  -)  -uexe(v+  -)2  +ueye(v+  -)3  v=0,l,2,...v0  (51) 

where  h  is  the  Planck’s  constant  and  c  is  the  speed  of  light.  Here,  denotes  ground  state 
vibrational  energy,  e2  denotes  first  excited  state,  and  so  on.  The  spectroscopic  constants 
are  given  by,56  cue=2358.57  cm-1,  uexe= 14.324  cm-1,  and  u>eye=-0. 00226  cm-1.  When 
u0=47,  the  value  of  energy  exceeds  the  N2  dissociation  energy,  9.86  eV. 

3.5  Vibrational  Bias  and  Depletion  in  Dissociation 

With  an  increase  in  flow  temperatures,  an  accumulation  of  internal  energy  leads  to 
dissociation  of  the  diatomic  molecules.  The  dissociation  process  in  the  continuum 


31 


(52) 


formulation  is  modeled  in  the  form  of  an  Arrhenius  equation 

k,(Te ff)  =  CVr>ffexp(fVTeff) 

where  the  rate  controlling  temperature,  Tcq  requires  an  appropriate  definition.  Present  day 
hypersonic  codes  have  often  used  the  Park  model  to  define  the  rate  controlling  dissociation 
temperature,  due  to  its  ease  of  implementation.  However,  the  model  has  been  shown  to  be 
inaccurate  for  a  wide  range  of  temperatures.  In  order  to  highlight  the  kinetic  implications 
of  nonequilibrium  relaxation,  dissociation  models  differing  in  the  vibrational  bias,  suitable 
for  a  wide  range  of  degrees  of  nonequilibrium  are  discussed.  The  present  study  examines 
the  behavior  of  these  continuum  dissociation  models  in  oxygen  and  nitrogen  baths 
maintained  at  a  gas  temperature  of  6,000  K  and  10,000  K,  respectively  for  the  vibrational 
temperature  ranging  between  2,000  and  12,500  K.  Brief  description  of  the  dissociation 
models  follows. 

The  Park  model57  expresses  the  dissociation  coefficient  in  terms  of  an  effective  temperature: 

Teff  =  (53) 

kdef[  =  kdeq(TeS)  (54) 


The  ladder  and  weak  bias  model  were  developed  considering  the  energy  in  the  vibrational 
quantum  states  of  the  dissociating  diatomic  molecule.  Dissociation  in  the  ladder  model  is 
assumed  to  take  place  only  from  the  last  bound  state  whereas  for  the  weak  bias  model  all 
quantum  states  contribute.  The  dissociation  rate  of  the  Ladder  model  can  be  written  as, 


n*vkv »_ 


> Continuum 


(55) 


where,  n*  is  the  vibrational  population  in  the  last  state  and  kv* -^continuum  is  the 
dissociation  rate  from  the  last  vibrational  state. 

The  weak  bias  model,  in  contrast  to  the  previous  model,  assumes  dissociation  to  proceed 

from  all  the  states,  viz.  /((i/w) — >Continuumi  Tl\k\ — + Continuums  ^2 kL2 — » Continuum  and  SO  Oil  to 
kv*  — > Continuum  • 

In  order  to  facilitate  a  comparison  of  the  global  rate  of  dissociation  under  different 
conditions  of  vibrational  bias,  it  is  useful  to  define  an  effective  dissociation  rate, 
kd^[T,  Tv],  where: 

V* 

kdeg[T,  Tv]Ntot  =  kdvnv  (56) 

v=0 

and  kdv  is  the  dissociation  coefficient  from  vibrational  state  v  and  is  assumed  to  be  a 
function  of  T  only.  By  expressing  the  functional  dependence  of  the  population  in  each 
vibrational  state  in  terms  of  a  vibrational  temperature: 


nv 


Ntot  exp 


(  —€v 


) 


Jvib 


(57) 


where  Zvn,  is  the  vibrational  partition  function,  ev  is  the  vibrational  energy  in  the  quantum 
level  v,  k  is  the  Boltzmann  constant,  Ntot  is  the  total  number  density,  and  Tv  is  the 
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vibrational  temperature,  one  can  express  the  kd  g.[ T,TV ]  for  the  ladder  and  weak  bias 
model  in  terms  of  the  equilibrium  coefficient  kd[T] .  In  equilibrium, 

few  PI  =  jr  =  f;  ^exp[-£”/*T1 


v=0  ^tot  „=0 


ZVib[T] 


For  the  Ladder  model, 

kd  JT,TV]  =  kdvt^-  =  kdeq[T]^L  =  kdeq[T }  exp[e,A  -  ^-)] 

I'tot  Tl yeq  1  1V  Zjyibyly  J 

A  simple  weak  bias  model  may  be  constructed  by  assuming: 

fefc  =  fejp]  exp[-(D  -  e„)/T] 


(58) 

(59) 

(60) 


where,  the  value  of  k0  is  established  by  requiring  consistency  with  the  equilibrium  rate 
coefficient: 

fcopl  =  (61) 

Substitution  for  kdv  yields  for  this  weak  bias  form: 


^deff 


kdeq[T]  Zvib[T] 
( v *  +  1)  Zvib[Tv] 


v  1 

Ntot  y  exp  [ev  ( — 

D=0  1 


(62) 


Note  that  nonequilibrium  values  of  the  effective  dissociation  coefficient  for  the  ladder  and 
weak  bias  forms  differ  substantially.  The  state-specific  dissociation  rates  compete  with  the 
V-V  and  V-T  rates  and  the  state  populations  are  depicted  in  the  excitation  phase. 

The  above-mentioned  Ladder  model  does  not  account  for  depletion  of  the  vibrational 
population.  In  the  present  study,  the  Ladder  model  with  the  effect  of  population  depiction 
was  also  implemented.  The  vibrational  population  depletion  on  account  of  energy  loss  of 
dissociation  was  modeled  according  to  the  approach  outlined  by  Osipov  and 
Stupochenko.36  This  model  was  developed  by  Josyula  and  Bailey  in  Refs.  38  and  40,  the 
study  accounted  for  effects  of  V-T  and  V-V  transition  rates.  An  analytic  expression  of  the 
resulting  reduction  in  the  equilibrium  dissociation  rate  was  derived  by  defining  (pv  to  give 
the  deviation  of  the  quasi-steady  distribution  from  an  equilibrium  Boltzmann  distribution. 
The  mass  density  in  vibrational  state  v  is: 

Pv(t)  =p(v)(l  +  <pv)  (63) 


The  deviation  of  the  population  of  level  v  from  equilibrium  is  given  by  ipv  which  relates  pv 
(the  state  mass  density  in  level  v)  to  (the  state  mass  density  in  level  v  at  equilibrium). 

It  is  seen  from  Eq.  63  that  when  ipv  <  0,  the  level  population  decreases  relative  to  the 
equilibrium  population. 

The  coupled  conservation  equations  used  the  Roe  approximate  Riemann  solver, 
implemented  in  finite  volume  formulation  by  computing  the  cell  interface  flux  as  a 
summation  of  wave  speeds,  similar  to  Eq.  38.  The  second  order  spatial  accuracy  was 
obtained  by  employing  the  MLISCL  approach  in  conjunction  with  the  minmod  limiter  to 
reduce  the  solution  to  first  order  accuracy  in  the  vicinity  of  strong  shock  waves,  as 
described  in  the  work  of  Josyula,  et  al.54  An  explicit  predictor-corrector  method  is  used  to 
advance  the  solution  in  time. 
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Figure  17:  Typical  Computational  Grid  for  Blunt  Body  Flow 

3.6  Results  and  Discussion 

Results  are  presented  in  three  sections.  In  the  first  section,  the  multitranslational 
temperature  comparisons  with  DSMC  solutions  of  the  shock  structure  in  translational  and 
rotational  nonequilibrium  for  Mach  1.2  to  5.0  shock  waves  will  be  presented.  This  will  be 
followed  by  the  flow  simulations  past  two-dimensional  blunt  bodies  for  Mach  numbers  from 
6.5  to  19.83  to  assess  the  effect  of  multiquantum  V-T  energy  exchanges  on  the  relaxation 
process  and  the  effect  of  vibrational  bias  in  the  dissociation  process.  A  typical 
computational  grid  used  for  the  solution  of  the  axisymmetric  blunt  body  flows  in  thermal 
and  chemical  nonequilibrium  is  shown  in  Figure  17. 

3.7  Multitranslational  and  Rotational  Temperatures  in  Shock 
Structure 

To  determine  the  predictive  accuracy  of  the  multitranslational  temperature  model,  the 
internal  structure  of  a  shock  wave  was  simulated.  Since  there  are  no  experimental  data  for 
the  two-translational  temperature  measurements,  the  present  calculation  will  be  compared 
with  the  DSMC  solutions.  Currently,  the  DSMC  is  the  most  reliable  method  to  obtain 
solutions  of  a  shock  structure  that  is  in  thermodynamic  nonequilibrium.  The  present 
calculations  use  a  monatomic  gas  with  the  atomic  weight  of  argon  and  a  Maxwellian 
(inverse  fourth  power)  inter-atomic  potential  (/ire/  =  2.515  x  10 ~5kg/ms,  Tref  =  273 K,  and 
uj  =  1).  All  cases  use  an  equal  space  mesh  of  200  grid  points.  Candler  et.  al.51  conducted  a 
study  of  shock  structure  with  multiple  translational  temperatures  using  a  continuum 
method  and  compared  with  DSMC  results.  The  continuum  study  of  Ref.  51  predicted  a 
much  thinner  shock  width  than  that  of  DSMC  predictions.  The  solutions  of  the  present 
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Figure  18:  Computed  Normalized  Density  and  Temperatures  for  a  Mach  1.2  Maxwellian 

Gas  Shock  Wave 


study  will  be  compared  with  the  DSMC  shock  structure  predictions  presented  in  Ref.  51. 

Figures  18  through  20  plot  the  Maxwellian  gas  shock  waves  for  free-stream  Mach  numbers 
between  1.2  to  5.0. 

The  normalized  temperature,  Tn  =  (T  —  T, ) / (T2  —  7\),  and  normalized  density 
Pn  =  (p  —  pi)/ (p‘2  —  pi )  are  plotted  versus  x/\\.  The  upstream  mean-free  paths  Ai  is 
dehned  as  Ai  =  2/x/piCi,  where  C\  is  the  mean  atomic  speed  at  the  upstream  conditions. 
The  temperature  is  resolved  in  the  components  perpendicular  (x  direction)  and  parallel  (y 
direction)  to  the  shock  wave.  From  these  figures,  we  can  see  the  excellent  match  between 
the  current  results  and  DSMC  solutions  for  all  Mach  number  cases  from  the  continuum 
(M  =  1.2)  to  the  highly  nonequilibrium  one  (M  =  5).  For  the  Mach  5  flow,  there  is  slight 
disagreement  with  the  DSMC  solution  in  the  upsstream  location.  The  density  and 
rotational  temperature  have  also  been  compared  with  the  experimental  data  of  Robben 
and  Talbot.58  The  computation  considered  a  pre-shock  temperature  of  9.15  K.  The 
measured  rotational  temperature  and  density  have  been  observed  to  be  accurately 
predicted  by  the  present  continuum  model. 

3.8  Effect  of  MultiQuantum  V-T  Exchanges  on  Vibrational 
Temperatures  in  Blunt  Body  Flows 

The  flow  conditions  for  the  results  shown  in  the  Figures  21  and  22  are:  M00=6.5,  Poo—^0 
Pa  and  Too=300  K,  for  pure  nitrogen  flow  past  an  infinite  cylinder  of  1  meter  radius. 
Figure  21  shows  temperature  distributions  in  the  flowheld  for  computations  assuming 
single,  double  and  six-quantum  V-T  energy  exchanges,  with  and  without  single  quantum 
V-V  energy  exchanges;  Figure  21a  along  the  stagnation  streamline  and  Figure  21b  along 
the  surface.  An  anharmonic  oscillator  with  40  quantum  levels  for  the  nitrogen  molecule 
was  assumed  and  the  thermal  nonequilibrium  contribution  was  coupled  to  the  fluid 
dynamics  equations.  The  results  indicate  that  the  influence  of  multiquantum  V-T 
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Figure  19:  Computed  Normalized  Density  and  Temperatures  for  a  Mach  2.0  Maxwellian 

Gas  Shock  Wave 


Figure  20:  Computed  Normalized  Density  and  Temperatures  for  a  Mach  5.0  Maxwellian 

Gas  Shock  Wave 
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Figure  21:  Comparison  of  Temperatures  Along  Stagnation  Streamline  and  Along  Surface 

for  Nitrogen  Flow  Past  Blunt  Body 


Single  Quantum  V-T 


Figure  22:  Population  Distribution  at  Select  Locations  Along  Stagnation  Streamline  for 

Nitrogen  Flow  Past  Blunt  Body 
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Figure  23:  Comparison  of  Vibrational  Biases  on  the  Effective  Dissociation  Rate  Coefficient 

transitions  on  the  translational  and  first-level  vibrational  temperatures  along  the 
stagnation  streamline  and  surface  are  negligible.  The  results  show  a  small  influence  of  the 
inclusion  of  single  quantum  V-V  exchanges  on  the  translational  temperature  along  the 
stagnation  streamline,  V-V  energy  transfers  speed  up  the  relaxation  process.  There  is, 
however,  a  more  pronounced  effect  of  V-V  exchanges  on  the  first  level  vibrational 
temperature.  Figure  22  shows  the  population  distribution  for  nitrogen  along  the  stagnation 
streamline  at  the  freestream,  pre-shock,  post-shock  and  stagnation  point  locations.  The 
negligible  influence  of  multiquantum  transitions  are  also  seen  in  these  results.  At  the  lower 
Mach  number  of  6.5  considered  in  the  present  study,  it  may  be  concluded  that 
multiquantum  V-T  exchanges  have  a  minor  effect  on  the  vibrational  relaxation  process. 

3.9  Effect  of  Vibrational  Bias  in  Blunt  Body  Flows 

Figure  23  shows  the  calculated  two-temperature  dissociation  rate  constant  kd  ^(T,TV)  for 
the  cases  of  oxygen  and  nitrogen  baths  maintained  at  T  —  6,  000  K  and  T  =  10,  000  K, 
respectively.  The  vibrational  temperature  varies  from  500  K  to  12,500  K.  This  represents 
the  basic  condition  that  arises  in  shock  layer  compression  flows,  where  the  internal  mode 
temperature  lags  that  of  translational  temperature.  The  comparison  of  different  biases  of 
the  dissociation  mechanism  reveals  large  differences  at  vibrational  temperatures  far  from 
equilibrium.  From  previous  studies  (See  for  example  Ref.  41)  it  is  conjectured  that  for  for  a 
Tv/T  less  than  0.5,  the  dissociation  occurs  from  the  lower  levels  and  for  Tv  ~  T, 
dissociation  can  occur  from  all  the  levels.  The  figure  shows,  for  T  >  Tv,  the  rate  coefficients 
of  the  Ladder  model  (signihying  strong  bias)  are  several  orders  of  magnitude  lower  than 
those  of  the  weak  bias  model.  For  T  <  Tv,  the  trend  is  reversed  and  the  Ladder  model  has 
a  higher  rate  coefficient  than  the  weak  bias  model.  The  results  from  the  Park  model  are 
shown  for  reference.  It  may  be  concluded  that  the  stronger  the  bias,  the  more  rapid  is  the 
decrease  in  the  rate  coefficient  for  Tv  <  T . 

Figure  24  shows  the  inverse  of  the  term  in  Equation  (63),  1  +  (pi  which  gives  the 
dissociation  reduction  factor  from  Park’s  equilibrium  rates  plotted  as  a  function  of 
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Dissociation  reduction  factor 

1  /  (1  +9,) 


Figure  24:  Dissociation  Reduction  Factor  for  Oxygen  and  Nitrogen  Due  to  Depletion 

temperature.  Note  that  when  <pi  <  0,  the  effect  is  to  deplete  the  population  in  the 
vibrational  states.  The  reduction  factor,  is  the  ratio  of  the  level  population 
distribution  at  equilibrium  to  the  nonequilibrium  state.  The  reduction  factor  which 
accounts  for  the  deviation  from  the  quasi-steady  state  distribution  can  range  from  1  to  3 
orders  of  magnitude  between  the  temperatures  of  5,000  K  to  10,000  K  at  which  there  is 
considerable  nitrogen  dissociation.  This  factor  was  applied  as  a  correction  to  the  Park’s 
dissociation  rates  to  account  for  the  depletion  effects  in  the  vibrational  levels,  ft  is 
important  to  note  that  the  reduction  factor  is  the  highest  below  2,000  K  but  there  is 
negligible  nitrogen  dissociation  at  this  temperature.  The  reduction  factor  again  peaks  at 
about  6,000  K  and  stays  high  till  12,000  K,  a  region  of  significant  nitrogen  dissociation  and 
at  higher  temperatures  begins  to  fall  and  reaches  a  factor  of  1  at  the  highest  temperature 
of  15,000  K.  The  dissociation  reduction  factor  in  oxygen  dissociation  also  shown  in 
Figure  24  up  to  a  temperature  of  10,000  K  shows  a  similar  trend  as  for  nitrogen. 

Figure  25  shows  the  effect  of  the  various  dissociation  models  including  the  ladder  model 
with  depiction  effects  on  stagnation  streamline  profiles  for  a  Mach  19.83  nitrogen  flow  past 
a  blunt  body.  The  Park  dissociation  model  captures  qualitatively,  [Figure  25(a)],  the 
expected  nonequilibrium  between  the  translational  and  vibrational  modes.  Shown  here  are 
the  temperature  profiles  with  depletion  effects;  the  large  depletion  (see  Figure  24  for 
dissociation  reduction  factor),  reduces  the  effective  dissociation  rate,  thereby,  increasing 
the  shock- standoff  distance.  Figure  25(b)  compares  profiles  of  atomic  nitrogen  production 
for  the  various  models.  Note  the  extended  region  of  near-equilibrium  conditions  prevailing 
in  the  shock  layer  for  this  strong  shock  case.  The  Ladder  model  with  depletion  differs  from 
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Figure  25:  Effect  of  Vibrational  Bias  and  Depletion  in  Blunt  Body  Flow  on  Temperature 

and  Dissociation  Rate) 


the  other  models,  though  all  predict  similar  values  for  the  peak  production  rate  in  the 
shock  layer  due  to  the  near  equilibrium  conditions.  It  is  concluded  that  the  Ladder  model 
(without  depletion)  and  weak  bias  models  give  an  lower  and  upper  bound  on  the 
dissociation  rate  in  the  narrow  region  of  shock  jump  location. 

At  a  lower  Mach  number  of  11.18  for  air  flow  past  a  blunt  body,  there  is  considerable 
oxygen  dissociation  with  negligible  nitrogen  dissociation.  Temperature  profiles  along  the 
stagnation  streamline,  Figure  26,  shows  the  relatively  small  influence  of  bias  and  the  large 
influence  of  vibrational  population  depiction  in  oxygen  dissociation,  (See  Figure  24  for 
dissociation  reduction  factor).  The  shock  standoff  distance  predicted  with  the  effect  of 
depletion  in  the  Ladder  model  exhibits  a  better  agreement  with  experiment,  Figure  27  than 
the  weak  bias,  Park,  and  ladder  models. 
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Figure  26:  Comparison  of  Translational  and  Vibrational  Temperature  Along  the 

Stagnation  Streamline 


Figure  27:  Comparison  of  Shock-Standoff  Distances 
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4  Conclusions 


A  high-order  finite-difference  method  has  been  developed  for  use  in  3-D  multidisciplinary 
applications.  The  key  ingredients  in  the  approach  are  Pade-type  differencing  and  filtering 
formulas.  The  method  has  been  extended  to  treat  the  complications  associated  with 
nonuniform  meshes,  nonlinearities  and  truncated  boundaries.  In  particular,  issues  related 
to  boundary  condition  implementation  and  instability  due  to  mesh  stretching  have  been 
overcome.  Special  care  is  exerted  to  develop  the  proper  metric  evaluation  procedures  to 
retain  the  advantages  of  higher  order  accuracy  in  curvilinear  and  dynamically  deforming 
3-D  meshes.  An  effective  overlap-based  strategy  is  combined  with  boundary  formulations 
to  facilitate  domain-decomposition  in  parallel  processing  strategies  for  arbitrary  meshes.  A 
wide  range  of  numerical  experiments  demonstrate  that  the  method  overcomes  the 
traditional  problems  which  had  limited  the  utilization  of  high-order  methods  to 
multidisciplinary  phenomena  in  practical  situations. 

For  high-temperature  physics,  a  first-principles-based  method  has  been  developed  to  solve 
the  interaction  between  internal  modes.  Numerical  simulations  have  been  completed  of 
steady  state,  hypersonic  flows  for  the  prediction  of  multitranslational  and  rotational 
temperatures,  and  vibrational  energy  distributions.  Separately,  a  new  gas  kinetic, 
multitranslational  temperature  model  for  a  monoatomic  gas  has  been  developed  by  solving 
equations  derived  from  the  Boltzmann  equation  with  a  first  order  Chapman-Enskog 
expansion  of  an  anisotropic  velocity  distribution  function.  The  model,  extended  to 
rotational  nonequilibrium,  has  been  employed  for  highly  nonequilibrium  shock  structures 
and  compared  with  DSMC  solutions  and  experimental  data.  The  nonequilibrium 
vibrational  energy  distributions  have  been  modeled  by  the  master  equation  and  the 
population  distributions  in  the  quantum  energy  states  of  the  diatomic  molecule  evaluated 
under  multiquantum  vibrational-translational  energy  exchanges.  The  dissociation  process 
considered  the  effect  of  vibrational  bias  and  depletion  on  the  equilibrium  dissociation  rate. 
Results  of  hypersonic  blunt  body  flow  have  compared  with  experimental  shock  stand-off 
distance.  It  has  been  shown  that  the  Ladder  model  with  population  depletion  provided  the 
most  favorable  comparison  with  experiment  compared  with  the  Park  model,  weak  bias 
model,  and  Ladder  model  without  depletion.  These  results  guide  a  better  understanding 
and  characterization  of  the  nonequilibrium  processes  prevalent  in  hypersonic  aerospace 
vehicle  flowfields. 
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a,  6,  C 
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Cp 
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Fvi  GVJ  H % 
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kdv 
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Kn 
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Tv 

U,  V,  w 
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V-T 

V-V 

x,y,z 

Zyib 

a,  (3 


coefficients  of  explicit  terms  in  the  compact  difference  formula 
coefficients  of  explicit  terms  in  the  compact  filter  formula 
specific  heat  at  constant  pressure 
total  energy  per  unit  mass 
total  specific  energy 

nondimensional  components  of  the  electric  field  vector 

molecular  velocity  distribution  function 

inviscid  vector  fluxes 

viscous  vector  fluxes 

anisotropic  distribution  function 

grid  point  index  number 

Jacobian  of  the  coordinate  transformation 

Boltzmann  constant 

dissociation  rate  from  state  v  [cm3/s\ 

Equilibrium  constant 

Knudsen  number 

Molecular  weight  of  species  i 

Mach  number 

molecular  weight,  kg /mole 

number  density 

Noncommunicating  overlap 

pressure,  N/m 2 

Prandtl  number,  0.73  for  air 

heat  flux  vector,  W 

vector  of  dependent  variables 

components  of  the  heat  flux  vector 

gas  constant 

reference  Reynolds  number,  PoowooV/hx> 
source  vector 

equilibrium  trans-rotational  temperature 
nondimensional  time 
translational  temperature 
vibrational  temperature 

nondimensional  Cartesian  velocity  components 

vibrational  quantum  numbers 

Vibrational-  Tr  ansla  t  ional 

Vibrational-Vibrational 

nondimensional  Cartesian  coordinates 

vibrational  partition  function 

coefficients  in  difference  formula 
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Subscripts 

max 

min 


coefficients  in  filter  formula 
quantum  level  energy 
ratio  of  specific  heats 
thermal  conductivity  coefficient 
upstream  mean  free  path,  m 
viscosity,  kg/m-s 
general  scalar  variable 
characteristic  time  scale 
Navier-Stokes  stress  tensor,  N/m 2 
characteristic  temperature  of  dissociation 
characteristic  temperature  of  vibration 
depletion  factor 

total  density,  total  molecular  population 
nondimensional  body-fitted  computational  coordinates 
metric  coefficients  of  the  coordinate  transformation 


maximum  value 
minimum  value 
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